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The  nonlinear  boundary  value  problem  for  the  propaga- 
tion of  random  surface  gravity  waves  in  an  ocean  of  finite 
depth  is  solved  correct  to  second  order  in  perturbation 
parameter.   A  nonlinear  second  order  interaction  kernel  is 
obtained  which  results  in  nonlinear  spectral  corrections  to 
a  linear  Gaussian  sea  spectrum  at  frequencies  which  are  the 
sums  and  differences  of  the  interacting  frequencies  of  the 
linear  spectrum.   The  trivariance  function  is  shown  to  be  a 
closed  statistical  measure  of  the  second  order  nonlinearities 
and  the  coarse  measure  of  the  trivariance  function,  i.e.,  the 
skewness,  is  used  to  determine  the  magnitude  of  the  second 
order  nonlinear  contributions.   An  algorithm  is  presented  for 
simulating  a  time  sequence  of  nonlinear  random  surface  gravity 
waves  correct  to  second  order  by  employing  the  fast  Fourier 
transform. 

The  normalized  cumulative  probability  distributions  of 
linear  and  nonlinear  simulated  realizations  are  compared  with 


measured  hurricane-generated  realizations  recorded  during 
Hurricane  Carla  in  the  Gulf  of  Mexico  by  Wave  Force  Project 
II.   The  simulations  are  synthesized  from  a  smoothed  mea- 
sured spectrum  and  from  a  Bretschneider*  spectrum  having 
equal  variance  and  best  least-squares  fit  to  the  smoothed 
measured  spectrum  in  order  to  demonstrate  the  effect  of  spec- 
tral shape  and  phase  angles  in  random  simulations  in  addition 
to  establishing  the  applicability  of  the  Bretschneider  spec- 
trum for  design. 

Resultant  wave  forces  at  the  55  foot  dynamometer  eleva- 
tion are  predicted  from  both  the  linear  and  nonlinear  simu- 
lations synthesized  from  both  spectra.   Kinematic  fields  are 
computed  by  the  digital  linear  filter  technique  modified  by 
a  vertical  coordinate  stretching  function  and  are  used  in 
the  Morison  equation  with  Dean  and  Aagaard  drag  and  modified 
inertia  coefficients.   Normalized  cumulative  probability 
distributions  and  wave  force  spectra  of  the  simulated  forces 
are  compared  with  measured  wave  force  realizations  and  with 
wave  force  realizations  predicted  by  filtering  the  measured 
sea  surface  realization. 


XIV 


CHAPTER  1 
INTRODUCTION 

As  the  design  of  offshore  permanent  pile -supported  struc- 
tures moves  beyond  the  four  hundred  foot  bottom  contour  into 
regions  of  greater  depth  and  as  the  mathematical  models  which 
describe  the  compliant  nature  of  the  pile-soil  interaction 
become  more  sophisticated,  the  requirement  to  perform  dynamic 
analyses  of  these  structures  becomes  more  critical.   Foster 
[46],  Edge  and  Meyer  [44],  Borgman  [20],  Malhotra  and  Penzien 
[83,84],  Nath  and  Harleman  [92],  Plate  and  Nath  [104],  Selna 
and  Cho  [114],  Mansour  and  Millman  [85],     Berge  and  Penzien 
[12],  Muga  and  Wilson  [91],  inter   alios,    have  presented  models 
for  the  dynamic  response  of  permanent  pile-supported  struc- 
tures to  random  forces.   The  random  forces  used  in  these 
studies  were  either  a  superposition  of  linear  Fourier  compon- 
ents with  random  phase  angles  which  were  uniformly  distributed 
between  (.-■n,-n')    or  a  strictly  periodic  Stokian  wave  of  finite 
amplitude.   The  purpose  of  this  study  is  to  provide  a  non- 
linear random  time  series  realization  of  a  surface  gravity 
wave  spectrum  in  an  ocean  of  finite  depth  and  to  utilize  this 
realization  as  a  forcing  function  to  a  filter  to  obtain  the 
kinematics  required  to  predict  pressure  forces.   The  emphasis 
will  be  on  the  simulation  of  a  random  time  series  vice  the 
determination  of  the  invariant  statistics  of  random  realizations 


In  Section  1,  the  most  important  results  and  conclusions 
obtained  by  others  in  related  investigations  toward  obtaining 
second  order  perturbation  corrections  to  the  sea  surface  are 
briefly  reviewed.   In  Section  2,  the  ef*fect  on  the  resulting 
sea  surface  profile  of  the  linear  superpositioning  and  the 
nonlinear  interaction  of  two  collinear  waves  is   examined  by 
means  of  a  simple  model.   In  Section  3,  the  application  of 
the  nonlinear  sea  surface  simulation  to  obtain  a  random  pres- 
sure force  time  series  is  briefly  outlined. 

1.   Related  Previous  Studies 

Tick  has  developed  a  method  for  generating  a  realization 
of  a  nonlinear  sea  surface  for  deep  water  [123]  and  for  water 
of  finite  depth  [125]  by  formally  perturbing  the  energy  spec- 
trum.  For  a  Gaussian  initial  estimate  to  the  linear  boundary 
value  problem,  the  perturbed  energy  spectrum  was  shown  to  con- 
sist only  of  even  powers  of  the  perturbation  parameter.   The 
energy  spectrum  was  then  used  as  the  measure  of  the  nonlinear- 
ities.   In  order  to  graphically  display  the  superposition  of 
the  linear  energy  spectrum  with  the  quadratic  energy  spectrum, 
the  ordinate  values  of  the  quadratic  spectrum  had  to  be 
greatly  expanded  near  the  origin  (even  in  the  cases  of  finite 
depth)  in  order  that  the  measure  of  the  quadratic  perturba- 
tion of  the  energy  spectrum  could  be  observed. 

In  his  initial  study  (cf.  Kinsman  [72]),  Tick  treated 
finite  amplitude  waves  correct  to  second  perturbation  order 


in   deep  water   and  he   obtained   a   correction   to   the    free    sur- 
face   displacements    given  by   the    following  expression: 

z   =   0, |x|«»,t>0  (1.1) 

where  the  alphabetical  subscripts  denote  partial  differen- 
tiation with  respect  to  the  independent  variables  indicated 
by  the  suffix  and  the  numerical  prefix  denotes  the  perturba- 
tion order  of  the  dependent  variable.   The  last  of  these  three 
terms  represents  the  Maclaurin  series  expansion  of  the  varia- 
tion of  the  free  surface  about  the  still  water  level  due  to 
the  presence  of  the  wave.   The  second  term  in  brackets  repre- 
sents the  correction  due  to  the  local  kinetic  energy  of  the 
water  particle  motions.   The  first  term  represents  the  second 
order  perturbation  correction  to  the  velocity  potential.   The 
second  order  perturbation  correction  to  the  velocity  potential 
in  deep  water  was  given  by  Tick  in  the  following  integral 
expression: 

00 

2<I>(x,z,t)  =  2i//G(a),cj')  exp{  |  (  |  a)|  a)+ |  to'  |  co' )  |  •  -}  •  ^f  (oj' )  • -,^f  (w)  • 

•exp    i{6  (w,x,t3  }-exp    i{6  (to' ,x, t)  }dwd(jo'  (1.2) 

where   the   nonlinear   interaction  kernel,    g((jo,u)'),    is    given  by 

G(tu,to')    =   03(0)0^'     -     hlh'l) (1.3) 

(co+o)')^    -    I  (|a)|a)+|to'  |(jj')  I 


Kinsman  [72]  corrected  the  separation  constant  originally  pub- 
lished by  Tick  [123]  in  order  that  the  velocity  potential 
satisfy  exactly  the  equation  of  continuity.   This  correction 
has  been  incorporated  in  Eq.  (1.2).   The  second  order  pertur- 
bation correction  to  the  free  surface  computed  from  Eq.  (1.1) 
by  Tick  was  also  expressed  by  an  integral  equation 

1 
2n(x,t)  =  ^   //H(w,(o')-  f(x,a3,  t)-,f'(x,a)' ,t)da)daj'       (1.4) 

6    -00  -'■  -'• 

where  the  nonlinear  interaction  kernel  is  given  by 

ff(a),w')  =  2(a)+a)')  •G(co,w' )  +  (w)^  -  j    |a)w'|+  j  (ii^'  (1.5) 

One  of  the  most  important  results  obtained  by  Tick  was 
that  the  contribution  to  the  perturbed  energy  spectrum  at  a 
given  frequency  by  the  second  order  perturbation  correction 
to  the  sea  surface  was  the  result  of  the  convolution  of  all 
first  order  perturbation  energy  densities  whose  sums  and  dif- 
ferences in  frequencies  contribute  to  the  given  frequency  of 
the  quadratic  spectrum.   The  expression  computed  by  Tick  was 

1 
28(0))    =   ^     ■/"   if(a)-aj'  ,03')  •,S(a)-oa')-,S(u)' )dw'  (1.6) 

nn         g      °°  nn  nn 

This    important   result  will  be    discussed   in  more    detail    in 
Chapter  2. 

Tick    [125]    later  extended  his    results    to  water  of   finite 
depth   and  presented  the    following  nonlinear   interaction 
kernel : 


(|k|k+Ik' |k')    tanh[d2(|k|k+|k'  Ik-)]    -    (03+0)')^ 

(1.7) 

where    d     =    (h/g)    and  the   wave   numbers   k   =   k((jj)    and  k'    =   k'(aj) 

are   solutions    to   the    dispersion   relation 

0)^   =  k^   tanh[d^k^]  (1.8) 

The  sign  of  the  second  term  given  in  Eq.  (1.7)  has  been  cor- 
rected from  the  typographical  error  in  Reference  [125].   Tick 
resolved  the  difficulties  in  obtaining  an  identity  between 
the  interaction  kernels  given  by  Eq.  (1.5)  and  Eq.  (1.7)  as 
the  water  depth  approaches  infinity  (i.e.,  h->+°°)  in  terms  of 
the  discontinuities  in  hypersurfaces  created  by  Dirac  delta 
functions  and  by  the  nonuniqueness  of  perturbation  expansions. 

Tick  [125]  demonstrated  a  comparison  between  a  linear 
realization  computed  from  the  Neumann- type  spectrum  [72]  with 
a  second  perturbation  order  realization  from  the  same  spectrum 
for  the  case  of  an  ocean  of  infinite  depth.   No  realizations 
were  available  from  the  finite  depth  case  due  to  the  computa- 
tional complexities  involved.   The  Fast  Fourier  Transform 
(FFT)  algorithm  [37],  which  was  not  generally  available  at 
the  date  Tick  published  his  work,  has  greatly  facilitated  the 
applications  of  quadratic  sea  simulations  by  reducing  signifi- 
cantly the  computational  difficulties.   No  published  compari- 
sons of  random  quadratic  sea  simulations  with  measured  wave 
data  seem  to  be  available. 


In   an   extraordinary   three-part   series,    Hasselmann    [54, 
55,56]    set    forth   a  powerful   theory    for  the  nonlinear  energy 
transfer   in    a   spectrmn   of  random  waves.      Later  extensions 
[58,59,61]    included   the   probabilistic   structure   of   the    initial 
conditions    imposed   through   the   nonlinear   free   surface   boundary 
conditions    as  well    as    a   general    theory    [60]    for   the   wave -wave 
scattering  processes    in   the   oceanic  wave   guide.      Recently 
[76,120],    portions    of  this   powerful  wave-wave   scattering 
theory  have   been   applied  to  measured   data  with   excellent 
agreement.      The    incredible  number  of  geophysical   problems 
which   are    covered  by   this    general    theory   offer  many   opportun- 
ities   for   comparison  with  measured   data  now  becoming   available. 

In   one    of   these   works,    Hasselmann    [54]    observed   that   the 
quadratic   energy   spectrum  which  was    employed  by   Tick  was   not 
closed  to   the    order  of  the   perturbation  parameter   chosen. 
Hasselmann   included   in  his    development   the   mean   lagged  prod- 
uct  of  the    initial    Gaussian  estimates   with   the    solution   to 
the    third  order  perturbation   correction   to   the   sea   surface 
realization.      Therefore,    in   order   to   measure   nonlinearities 
by   a   formal   perturbation   of  the   energy   spectrum,    the   nonlinear 
contributions    to   the    sea   surface    realizations   must   be    com- 
puted  from  the  boundary   value   problem  correct    to   third   order 
in   the   perturbation  parameter   if   a  perturbation   of   the   energy 
spectrum   is    to  be   utilized   as    the   measure   of  the   nonlineari- 
ties.     The   bookkeeping  problems   of  perturbation  expansion   orders 
are    discussed  by  Hasselmann   et_  al.     [57]    and  equations   which 
aid   in   this   bookkeeping   are   given.      Briefly,    the   problem   is 


if  the   energy   spectriom    (or,    equivalent ly,    the   autocovariance 
ftinction)    is    formally  perturbed  in   a  power  series   expansion 
in  the   following   forms: 

CO  * 

E(w)    =   Z   ^E(a))  (1.9a) 

n 
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y(t)   =  Eii  ^nCt+T)-E  ^n*Ct)}  (1.9b) 
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the   relationship  between   the   energy  perturbation   orders    and  the 

time   series   perturbation   orders   must  be    complete    and   closed. 

An   intermediate    result   obtained  by  Hasselmann    [54]    in 
solving   the   problem  of  determining   the    rate   of  energy   transfer 
between   discrete    components    in   random  gravity  wave      trains    in 
water  of   finite   depth   included   a  nonlinear   interaction  matrix 
which  may   also  be   used  to   construct   nonlinear   seas.       In  his 
power   series    expansion   of   the   energy   spectriim   for   the    sea   sur- 
face,   Hasselmann   included   all  perturbation   orders    omitted   in 
the   expansion   used  by  Tick. 

These   studies   by   Tick   and  Hasselmann  have    introduced 
methods    for   simulating   random  nonlinear   realizations  of   surface 
gravity  waves.      Other   investigations    related  to   nonlinear   sea 
simulations   have   been   directed  toward  describing   the   effects    of 
nonlinearities    on   the    statistics    of  the    distributions    of  the 
realizations.      In  one   of  these   statistical    investigations, 
Longuet-Higgins    [79]    showed   that   the   Gram-Chalier   distribution 
more    closely    fits   measured  sea   surface    recordings    and,   more- 
over,   the   main   distinction  between   the    Gaussian   and  Gram- 
Chalier   distributions    are   terms   proportional    to   the   measure 
of  the   skewness.      Longuet-Higgins    [79]    also   obtained  the 


important  result  for  a  discrete  spectrum  of  waves  first  given 
by  Tick  [12  3]  for  a  continuous  spectrum  of  waves  that  the  quad- 
ratic energy  distribution  is  the  result  of  the  summation  of 
all  discrete  energy  components  whose  sums  and  differences  in 
phases  contribute  to  the  same  spectral  frequency. 

Cartwright  and  Longuet-Higgins  [32]  have  also  demonstrated 
the  effect  of  the  bandwidth  of  a  spectrum  on  the  probability 
distribution  of  nonlinear  seas  with  nonzero  skewness.   There- 
fore, from  both  the  statistical  (measured  moments)  and  probabil- 
istic (characteristic  and  probability  density  functions)  results, 
the  best  measure  of  the  nonlinearities  and  the  non-Gaussian 
distribution  of  a  sea  simulation  is  the  third  statistical 
moment.   Furthermore,  it  will  be  shown  that  the  boundary  value 
problem  need  only  be  solved  to  second  perturbation  order  in 
order  to  measure  the  lowest  order  contribution  to  both  the  non- 
linearities  and  non-Gaussian  distribution. 

The  omission  of  certain  perturbation  orders  in  the  power 
series  expansion  of  the  energy  measure  raises  a  closure  prob- 
lem similar  to  that  encountered  in  turbulence  (cf.  Batchelor 
[4],  Kampe  de  Feriet  [70],  or  Lumley  [81]).   Attempts  to  re- 
solve the  closure  problem  as  well  as  the  irreversability  para- 
dox of  the  initial  value  problem  resulted  in  a  series  of  papers 
by  Benney  and  Saffman  [7],  Hasselmann  [61],  Saffman  [113], 
Benney  [8],  Benney  and  Newell  [9,10,11],  Newell  [94],  Ablowitz 
and  Benney  [1],  Ablowitz  [2],  and  Birkhoff,  Bona  and  Kampe  de 
Feriet  [15]  which  eloquently  describe  the  physics  of  these  two 
problems  in  very  sophisticated  mathematical  terms.   Intermediate 


results  obtained  in  reaching  their  final  objectives  include 
nonlinear  interaction  matrices  (for  the  discrete  wave  spectrum 
interactions)  or  nonlinear  interaction  kernels  (for  the  con- 
tinuous wave  spectrum  interactions)  which  may  be  used  to  con- 
struct nonlinear  sea  simulations.   The  sophisticated  mathe- 
matics utilized  in  these  papers  contribute  to  the  subordina- 
tion of  the  practical  applications  of  their  intermediate  re- 
sults to  the  problem  of  simulating  nonlinear  seas. 

The  results  obtained  in  this  study  circumvent  the  closure 
problem  by  appealing  to  the  trivariance  function  as  a  measure 
of  second  order  nonlinearities  and  circumvent  the  irreversi- 
bility paradox  by  appealing  to  stationarity  and  ergodicity. 

Finally,  Borgman  [20]  has  developed  a  technique  which  is 
based  on  the  statistical  distributions  of  time  series.   The 
procedure  involves  the  transformation  of  a  time  series  which 
is  normally  distributed  to  a  time  series  which  is  gamma  dis- 
tributed with  zero  mean  by  a  set  of  polynomials.   While  this 
technique  is  extremely  attractive  for  its  computational  effi- 
ciency, it  is  not  based  on  the  hydrodynamic  equations  of  motion, 


2.   Linear  Superposition  and  Nonlinear 
Wave -Wave  Interactions 


Whitham  [130]  has  developed  an  extensive  number  of  appli- 
cations of  the  variational  technique  for  the  average  Lagran- 
gian  which  demonstrate  the  wave-wave  interactions  of  nonlinear 
dispersive  waves.   The  important  effects  of  both  amplitude  and 
frequency  dispersion  are  well  covered.   However,  a  simple 
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example  of  the  separate  effects  of  linear  superposition  and  of 
nonlinear  interaction  for  two  dispersive  wave  components  which 
is  relatively  free  of  mathematical  complexities  is  instructive 
in  distinguishing  the  relationships  between  the  two  effects. 

Formulation  of  the  complete  boundary  value  problem  for  a 
continuous  spectrum  of  dispersive  surface  gravity  waves  will  be 
deferred  until  Chapter  2.   However,  the  results  of  the  formal 
solution  in  deep  water  may  be  summarized  to  provide  a  simple 
illustration  of  both  the  linear  superposition  and  the  nonlinear 
interaction  of  a  pair  of  discrete  wave  components  from  a  spec- 
trum which  are  collinear  and  propagate  in  the  positive  x  direc- 
tion.  The  resulting  sea  surface  profile  may  be  represented  by 

n(x,t)  =  -^C(x,t)  +  2S(^.t)  (2.1) 

where  -■  C  is  the  linear  first  order  contribution  and  j^   i^  the 
nonlinear  second  order  contribution.   For  a  spectrum  of  two 
collinear  wave  components,  the  first  order  contribution,  ,?, 
represents  the  linear  superposition  of  the  two  waves  and  the 
second  order  contribution,  j^y    represents  the  contribution  of 
nonlinear  products  of  the  two  waves.   The  linear  superposi- 
tioning  of  the  first  order  contributions  will  be  examined 
briefly  first. 

2. a.   Linear  Superposition 

Consider  a  simple  line  spectrum  consisting  of  only  two 
spectral  components  which  are  separated  by  a  small  frequency 
differential,  2Aa,  and  by  a  small  wave  number  differential, 
2A^.   Although  adjacent  spectral  components  usually  differ  in 
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the  magnitude  of  their  amplitudes  (i.e.,  the  spectrum  is  not 
flat),  the  two  adjacent  components  will  be  considered  to  have 
the  same  amplitude  equal  to  H/2,  in  order  to  simplify  this 
illustration.  * 

The  deep  water  wave  profile  resulting  from  the  linear 
superposition  of  these  two  waves  is  given  by 

1^  =  ^+  -^  11-  (2.2) 

where  the  two  components  are  given  by 

T]^   =  J   cos{a^-t-ic^-x-i|;^}  =  ^  cos  e^  (2.3a) 

n_  =  -J  cos{a_-t-1c_-x-i|;_}  =  ^  cos  e_  (2.3b) 

where 

a^.  =  a  ±  Aa  (2.4a) 

ic"^  =  ic  ±  Ale  (2.4b) 

^+    =    ^   ±   A^  (2.4c) 

Substituting  Eqs.  (2.3a,b)  into  Eq.  (2.2)  and  combining  gives 
the  following  familiar  profile: 

■j^C  =  H  cos{(Aa)t-(A^)-x-Ai/j}  cos{at-^-x-i|;}  (2.5) 

This    resulting  wave    contains    an   amplitude  which   is   modulated 
by   the    differences  between   the   phases    of  the    two   components 
and  which  propagates    at   the   group   velocity,    Aa/|AicI.      The   pro- 
file  of   the    resulting  wave   propagates    at   the  wave    speed,    a/|ic| 
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The  resulting  wave  is,  of  course,  the  well  known  beat  profile 
which  is  the  result  of  the  superposition  of  two  linear  dis- 
persive waves. 

There  are  two  time  and  space  scale's  involved  in  the  beat 
profile  given  by  Eq.  (2.5).   One  is  the  scale  of  the  envelope 
modulation  which  is  spatially  long  and  temporally  slow  com- 
pared to  the  scale  of  the  profile.   The  relatively  longer  and 
slower  scales  are  related  to  periodicities  of  Ak  and  Aa , 
respectively,  while  the  relatively  shorter  and  faster  scales 
of  periodicities  are  related  to  k  and  a,  respectively. 

2.b.   Nonlinear  Interactions 


The  formal  derivation  of  the  second  order  contributions 
in  Eq.  (2.1)  will  be  given  in  Chapter  2.   However,  for  the 
purposes  of  this  illustration,  consider  the  following  two 
linear  velocity  potentials  for  the  deep  water  wave  components 
given  by  Eqs .  (2.3a,b): 

<D^  =  ^  exp{k^-z}  sin{a^'t-i<^-x-ij;^}  =  ^  exp{k^-z}sine^  (2.6a) 

$_  =  ^  exp{k_-z}  sin{ajt-ic_-x-ii;_}  =  ^  exp{k_-z}sine_  (2.6b) 

with   the    linear   dispersion   equation   for  deep  water  given  by 

toy   =   gk,  (2.7) 

Equation  (1.1)  indicates  that  the  second  order  contribu- 
tion to  the  wave  profile  is  a  function  of  the  temporal  deriva- 
tive of  the  second  order  velocity  potential.   Using  a  dis- 
crete form  of  the  second  order  velocity  potential  from 
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Chapter  2  and  substituting  the  appropriate  derivatives  o£ 
Eqs.  (2.6a,b)  into  Eq.  (1.1)  with  the  sign  of  the  first  term 
changed  from  a  minus  to  a  plus  results  in  the  following 
second  order  contributions:  * 

2  2  2 

2^  =  ^k^-cos  29^  +  ^k-cos  20_  +  ^k  cos  2  {at-Ic-x-^;} 

4a .  •  a 


S^  K-a.)2 


1 ^1 


lCa^-aJ^-g|k^-k_| 


cos  2{(Aa)t-(Ak) -x-Aif)} 
(2.8) 


The  first  two  terms  in  Eq.  (2.8)  represent  the  nondispersive 
Stokian  waves  resulting  from  the  self  interactions  of  the  two 
spectral  components  given  in  Eqs.  (2.3a,b).   These  two  Stokian 
waves  contribute  to  the  skewness  of  the  spectral  components 
by  adding  both  to  their  troughs  and  crests.   The  third  term 
in  Eq.  (2.8),  which  represents  the  result  of  the  cross  prod- 
uct of  the  interactions  between  the  two  spectral  components, 
propagates  at  the  profile  phase  speed  for  the  linear  wave 
given  by  Eq.  (2.5)  but  with  twice  the  phase.   This  wave  will 
contribute  to  the  skewness  of  the  resulting  linear  profile. 
The  final  term  in  Eq.  (2.8)  represents  the  result  of  the 
cross  product  of  the  interactions  between  the  two  spectral 
components  which  propagates  at  the  envelope  phase  speed, 
Aa/|Ak|,  but  with  twice  the  phase.   The  contribution  which 
this  wave  makes  to  the  envelope  depends  on  the  sign  of  the 
term  in  the  bracket.   For  a^>a_y    the  expression  becomes  the 
following: 

2 
^  aj-[l-(a_/a^)^]  cos  2{  (Aa)  t- (Alt)  •  x-A^;-  j}         >  (2.9) 
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This  wave  is  180°  out  of  phase  with  respect  to  the  envelope 
phase.   Its  maxima  occur  under  the  nodes  of  the  envelope  and 
its  minima  occur  under  the  antinodes  of  the  envelope.   This 
wave  also  has  a  nonzero  value  when  averaged  over  either  the 
faster  temporal  or  shorter  spatial  scales  of  the  resultant 
profile  phase  given  by  Eq.  (2.5).   Longuet-Higgins  and  Stewart 
[77]  discuss  numerous  effects  which  result  from  this  wave  com- 
ponent. 

It  can  be  shown  that  the  right  side  of  the  boundary  con- 
dition equation  which  determines  the  second  order  velocity 
potential  in  deep  water  reduces  to  the  temporal  derivative  of 
twice  the  kinetic  energy.   The  second  order  contribution  to 
the  resulting  wave  profile  given  in  Eq.  (2.8)  above  may  be 
seen  from  Eq.  (1.1)  to  be  proportional  to  the  temporal  deriva- 
tive of  this  second  order  velocity  potential  and  to  the  kine- 
tic energy.   Since  these  two  equations  demonstrate  that  the 
kinetic  energy  plays  such  an  important  role  in  the  second 
order  contributions,  it  is  of  value  to  examine  this  term  in 
some  detail. 

For  the  case  of  deep  water  waves,  the  right  side  of  the 
combined  free  surface  boundary  condition  reduces  to  the  tem- 
poral derivative  of  twice  the  kinetic  energy  which  may  be 
expressed  by 

2^tt  ^  2  2*z  =  It  '^1*1^  '   ^  =  °'  '^'^"'  ^-°  ^^'^^^ 

Substituting  the  appropriate  spatial  derivatives  of  Eqs. 
(2.6a,b)  into  the  right  side  yields 
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t  1^*1^  =  It 


5*+    2  '*+    2  "      2 


8x 


3$ 


+    2 


^ 


3x 


9^ 
9z 


9$_ 


;    z=0,    |x|<",    t>0 

(2.11) 
It   may  be    seen  by   inspection   of  Eqs.    (2.6a,b)    that   the   verti- 
cal  and  horizontal   velocity   components    are    90°    out-of-phase 
and,    consequently,    the    sum  o£   their   squares   evaluated   at    z=0 
is    a   constant  which   is    independent   of  time.      These    represent 
Stokian   self-interactions    and   the   temporal   derivative    of   the 
first   term  in  brackets    in   Eq.     (2.11)    vanishes    in   deep  water. 

Substituting   the    appropriate   spatial   derivatives    of  Eqs. 
(2.6a,b)    into   the   nonvanishing   temporal   derivative    in   Eq. 
(2.11)    yields 


9t 


nV 


k_^k_  [  (Aa)sin{0^-e_}    +   a-sin{  e^  +  e_  }  ]  + 
_+k^'k_[CAa)sin{e^-e_}    -   a-sin{e_^+e_ }  ] 


(2.12) 


For  a  collinear  spectrum,  it  is  easily  seen  that  only  the 
differences  in  phases  result  in  nonvanishing  contributions. 
Moreover,  for  Stokian  second  order  waves  (i.e.,  Aa=AkEO),  the 
right  side  vanishes  and  the  linear  velocity  potential  given 
by  the  sum  of  Eqs.  (2.6a,b)  is  exact  to  second  order. 

For  the  interaction  of  a  directional  spectrum  of  deep 
water  waves,  the  contributions  from  the  horizontal  velocity 
fields  are  decreased  by  a  multiplicative  constant  proportional 
to  cos{a^_}  where  a^_  is  the  azmuthal  difference  between  the 
directions  of  propagation  of  the  two  wave  components  given  by 
Eqs.  (2.3a,b).   The  multiplicative  constants  in  Eq.  (2.12)  no 
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longer  may  be  factored  and  the  contributions  from  the  sums  of 
the  phases  of  the  interacting  wave  components  no  longer  cancel 
Consequently,  a  finite  contribution  to  the  second  order 
velocity  potential  from  both  the  sums  and  differences  of  the 
phases  of  the  horizontal  and  vertical  velocity  fields  is 
obtained. 

The  important  differences  between  strictly  periodic 
Stokian  interactions,  collinear  spectral  interactions,  and 
directional  spectral  interactions  are  demonstrated  simply  for 
the  deep  water  case.   For  the  case  of  finite  depth,  additional 
algebraic  factors  are  encountered  which  detract  from  its  use- 
fulness as  a  simple  illustrative  example. 

3.   Simulating  Wave  Forces  by  the 
Digital  Linear  Filter  Technique 

One  of  the  practical  applications  for  the  simulation  of 
a  random  sea  surface  realization  is  to  obtain  a  random  force 
input  for  the  dynamic  analyses  of  permanent  pile-supported 
structures.   Reid  [106]  has  developed  a  method  for  the  linear 
filtering  of  a  sea  surface  realization  in  order  to  obtain  the 
kinematics  required  to  compute  pressure  forces  by  the  Morison 
equation.   Wheeler  [129]  used  this  technique  to  compute  pres- 
sure force  coefficients  from  measured  hurricane  realizations 
and  obtained  very  low  errors  between  the  measured  and  pre- 
dicted peak  pressure  forces  at  varying  elevations  along  an 
instrumented  platform  piling.   Another  application  [67]  of 
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this  technique  using  hurricane  realizations  yielded  mean 
square  errors  for  pressure  forces  computed  over  the  crest 
portion  of  the  realizations  which  compared  favorably  with  mean 
square  errors  from  pressure  forces  computed  using  the  Dean 
Stream  Function  [39].   Since  these  relatively  successful 
applications  of  the  digital  linear  filter  technique  to  analy- 
sis involved  measured  sea  surface  realizations  with  finite 
skewness,  the  extension  of  linear  digital  filtering  to  design 
should  involve  simulated  nonlinear  realizations  with  nonzero 
skewness. 

The  theory  for  nonlinear  sea  similations  will  be  devel- 
oped in  Chapter  2.   Selected  realizations  will  be  filtered 
in  Chapter  3  to  obtain  the  kinematics  required  to  compute 
pressure  forces.   Pressure  forces  computed  from  nonlinear 
realizations  by  the  linear  digital  filter  will  be  compared 
with  measured  pressure  forces  recorded  from  hurricane  gener- 
ated waves. 


CHAPTER  2 

THEORY  OF  RANDOM  NONLINEAR 
WAVE -WAVE  INTERACTIONS 


Section  1  of  this  chapter  briefly  outlines  the  proba- 
bility assumptions  to  be  employed  in  formulating  the  random 
sea  problem.   Section  2  introduces  the  notations  to  be  used 
in  describing  random  functions  and  gives,  without  proofs, 
some  of  the  mathematical  transforms  to  be  used  in  developing 
the  theory  for  random  surface  gravity  waves  correct  to  second 
order.   Section  3  defines  the  third  order  statistical  moment 
which  will  be  utilized  as  a  measure  of  the  nonlinearities  of 
the  random  sea  simulation  and  shows  that  this  statistical 
measure  is  closed  at  second  order  in  the  perturbation  param- 
eter.  Section  4  formulates  the  random  boundary  value  problem 
and  solves  for  the  first  and  second  order  random  measures  and 
spectral  response  functions  required  to  simulate  nonlinear 
random  sea  realizations.   Section  5  returns  to  the  nonlinear 
statistical  measure  of  the  second  order  correction  to  the  sea 
surface  realization  which  was  developed  in  Section  3  and  dis- 
cusses its  Fourier  transform.   Section  6  computes  the  second 
order  spectral  density  of  the  sea  surface  realization  and 
shows  that  it  is  a  convolution  of  all  first  order  spectral 
densities . 
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1.   Random  Differential  Equations 

Many  physical  processes  demonstrate  random  variations 
which  cannot  be  modeled  by  physical  laws  that  may  be  solved 
by  deterministic  mathematics.   Wiener  [133]  gives  several 
examples  with  applications.   These  random  variations  may  enter 
into  the  mathematical  equations  which  are  used  to  approximate 
the  random  physical  process  either  through  (1)  the  dependent 
variables  or  (2)  the  coefficients  in  the  differential  equa- 
tions or  (3)  through  both  methods.   In  addition,  the  random 
dependent  variable  may  be  a  random  function  of  still  another 
random  function  through  a  filtering  process  (cf.  Soong  [116], 
Jazwinski   [68],  or  Sveshnikov  [119]). 

In  any  event,  the  description  of  a  physical  process  as 
"random"  requires  an  estimate  of  the  probabilistic  structure 
of  the  random  function.   When  these  random  functions  are 
determined  by  differential  equations,  there  are  at  least  four 
convergence  questions  which  must  be  satisfied  for  each  single 
convergence  question  in  the  deterministic  differential  equa- 
tion (cf.  Soong  [116],  Jazwinski  [68],  or  Sveshnikov  [119]). 

One  of  the  most  well  known  probabilistic  structures  is 
the  Gaussian  density  and  a  Gaussian  probabilistic  field  is 
frequently  assumed  in  many  hydrodynamic  problems.   Chorin 
[34]  describes  many  of  the  implications  of  assuming  a  Gaus- 
sian field  for  random  flow  problems  in  hydrodynamics.   Rosen- 
blatt [111]  and  Kampe  de  Feriet  [71]  address  specifically  the 
probabilistic  questions  posed  by  the  Gaussian  field  assumption 
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for  the  random  linear  sea  problem.   The  probabilistic  struc- 
ture of  many  physical  problems  are  frequently  assumed  to  be 
known  and  the  assumption  of  a  Gaussian  field  for  the  random 
linear  sea  will  be  used  herein.   A  brie*f  discussion  of  the 
probabilistic  assumptions  used  is  unavoidable  in  order  that 
the  limitations  and  interpretations  of  the  probability  struc- 
ture of  the  applications  of  the  results  to  be  derived  may  be 
fully  understood. 

The  solution  to  the  random  sea  problem  is  a  random  func- 
tion which  is  determined  by  a  partial  differential  equation 
(viz.,    the  Bernoulli  equation)  in  which  one  of  the  dependent 
variables  is  the  velocity  potential  which  is  also  a  random 
function.   The  random  velocity  potential  is  determined  from 
a  boundary  value  problem  in  which  the  coefficients  of  the 
governing  field  equation  and  of  the  boundary  conditions  are 
deterministic.   The  random  velocity  potential  is  a  random 
function  of  a  random  variable  called  the  phase  angle.   The 
probability  density  function  of  the  random  phase  angle  is 
assumed  to  be  independently  and  identically  uniformly  dis- 
tributed between  {-tt,7t}.   Papoulis  [96,97]  demonstrates  that 
a  sum  of  random  functions  which  depends  on  an  identically 
imiformly  distributed  random  variable  tends  very  rapidly 
toward  a  Gaussian  distribution  even  for  cases  when  the  sum 
includes  relatively  few  random  variables.   The  rigorous  proof 
of  the  tendency  toward  a  Gaussian  distribution  by  random 
fields  as  the  number  of  the  random  variables  increases  with- 
out limit  may  be  found  in  Rosenblatt  [112]  by  means  of  the 
Central  Limit  Theorem. 
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The  three  fiondamental  convergence  theorems  required  in 
deterministic  calculus  involve  continuity,  differentiation 
and  integration.   In  random  differential  calculus,  each  of 
these  three  fundamental  theorems  must  be  shown  to  converge 
(1)  in  probability,  (2)  in  mean  square,  (3)  in  distribution, 
and  (4)  in  "almost  sure"  probability  (cf.  Soong  [116]).   Mean 
square  calculus  will  be  used  to  determine  the  convergence  of 
functions  and  mathematical  operators.   The  linear  solution  to 
the  random  sea  problem  will  be  assumed  to  be  second  order 
stationary  and  homogeneous  in  the  statistical  moments.   This 
assumption  of  stationarity  and  homogeneity  renders  the  proba- 
bility distribution  independent  of  time  and  space,  respectively 
(cf.  Yaglom  [135]).   The  assumption  that  a  random  process  is 
second  order  stationary  or  homogeneous  has  been  widely  studied 
and  exploited  in  terms  of  the  correlation  theory  (cf.  Jenkins 
and  Watts  [69]).   The  tendency  toward  a  Gaussian  distribution 
of  a  sum  of  random  functions  which  are  identically  and  uni- 
formly distributed  between  {-v,tt}    results  in  the  probability 
structure  of  the  random  linear  sea  process  being  completely 
described  by  its  mean  and  variance.   The  complexities  associ- 
ated with  demonstrating  "almost  sure"  probability  are  dis- 
cussed by  Bartlett  [3]  and  will  not  be  postulated  for  the 
random  sea  problem.   The  determination  of  the  evolution  of 
the  probability  density  function  for  the  random  nonlinear 
sea  will  also  not  be  attempted  (cf.  Jazwinski  [68]). 

The  complete  formulation  of  the  covariance  equations  for 
the  random  linear  sea  problem  using  correlation  theory  may  be 
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found  in  Snyder  [120].   The  distinctions  between  and  simplifi- 
cations of  the  solutions  to  the  linear  sea  problem  for  the 
probability  density  functions  using  probability  theory  versus 
the  covariance  functions  using  correlation  theory  are  demon- 
strated by  comparing  the  papers  by  Rosenblatt  [111]  and  Kampe 
de  Feriet  [71]  with  the  paper  by  Snyder  [120]. 

2.   Introduction  to  Random  Functions 

Kinsman  [72;  cf.  Chs.  7  and  8]  develops  a  general  expres- 
sion for  the  three-dimensional  energy  spectrum  for  a  random 
sea  surface  realization.   His  development  involves  the  three 
concepts  of  (1)  Fourier  analysis,  (2)  stochastic  processes  and 
probability,  and  (3)  the  hydrodynamics  of  the  random  sea.   In 
the  general  derivation  (Kinsman  [72],  Ch.  7),  each  of  the 
three  concepts  are  developed  separately  in  order  that  the 
effects  of  modifications  of  a  single  concept  may  be  examined 
only  in  terms  of  the  one  concept  which  it  modifies.   The 
general  model  is  then  explicitly  defined  in  terms  of  both  a 
deterministic  and  a  stochastic  model  in  which  the  interdepen- 
dencies  between  the  three  concepts  are  examined  (Kinsman  [72], 
Ch.  8). 

This  procedure  utilized  by  Kinsman  [72]  to  develop  a 
model  for  the  energy  spectrum  of  a  random  function  will  be 
followed  in  obtaining  a  solution  for  the  second  order  contri- 
bution to  a  random  sea  surface  realization.   However,  the  nota- 
tional  conventions  used  to  describe  random  functions  differ  in 
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some  respects  from  reference  to  reference.   In  order  to  avoid 
possible  confusion,  as  well  as  for  ease  of  reference,  the 
mathematical  expressions  and  definitions  to  be  used  in  devel- 
oping the  second  order  theory  of  a  random  nonlinear  sea  are 
given  in  this  section. 

The  general  definitions  for  random  functions  and  for  the 
correlation  theory  functions  with  both  time  and  space  argu- 
ments are  given  by  Kinsman  [72;  Chs .  7  and  8]  and  by  Phillips 
[100;  Ch.  4].   Definitions  for  functions  with  only  temporal 
arguments  are  given  by  Sveshnikov  [119]  or  Yaglom  [135].   The 
solution  developed  in  Section  4  of  this  chapter  is  initially 
a  solution  to  a  spatial  boundary  value  problem.   Therefore, 
the  definitions  given  below  are  for  functions  having  only 
spatial  arguments. 

A  random  function  which  is  a  continuous  function  of  time 
or  space  is  called  a  random  process.   A  random  function  which 
is  a  discrete  function  of  time  or  space  is  called  a  random 
sequence  (cf.  Sveshnikov  [119],  p.  11).   Any  random  process 
may  be  represented  by  the  following  Fourier-Stielt jes  integral: 

00 

C(x)  =  /  exp  i{i^-x}-d?(lt)  (2.1) 

-00 

where  C(k)  is  a  random  measure  with  orthogonal  increments. 

For  the  process  represented  by  Eq.  (2.1)  to  be  a  real  function, 

the  following  complex  conjugate  relationship  must  be  satisfied: 

dC(-it)  =  d^*(t)  (2.2) 
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A  good  discussion  of  the  distinctions  between  Fourier-Stieltjes 
integrals,  Stieltjes  integrals,  and  Fourier  integrals  may  be 
found  in  Yaglom  [135;  Ch.  2,  Sec.  9].   In  Eq.  (2.1)  and  in 
those  which  follow,  the  parameter  of  the  argument  is  space 
but  a  representation  for  an  argument  in  which  the  parameter 
is  time  may  be  similarly  defined.   The  autocovariance  func- 
tion for  a  process  with  random  measure  is  given  by 

Y(x,r)  =  E{c(x)-C(x+r)}  (2.3) 

For   a  process   which    is    spatially  homogeneous    (or   temporally 
stationary) ,    the    autocovariance    function   depends    only   on   the 
spatial    lag   r    (or   temporal    lag   t)    and  has    spectral    represen- 
tation  given  by   Khinchin's    Theorem    (cf.    Sveshnikov    [119],    Ch. 
II,    Sec.    10,    or  Yaglom    [135],    Ch.    2,    Sec.    10): 

00 

Y(r)    =/    exp    i{i^-r}    ds  (t)  (2.4) 

where   the    spectral    distribution    function,    s(k),    is    related  to 
the    random  measure,    ^ (k) ,    for   a  homogeneous   process   by 

£{dC(^)d?(ic')}    =   s(ic')-6(ic'+ic)d]^'dic  (2.5) 

where  6(-)  is  the  Dirac  delta  distribution  (Friedman  [49]). 
This  distribution  may  be  defined  in  the  space  domain  by 

-  oo 

(27r)^6(±x)    =/   exp    i{±ic-x}dic  (2.6) 

-  00 

and   in   the  wave   number   domain  by 

oo 

(2Ti)^6(±ic)    =/   exp   i{+ic-x}dx  (2.7) 
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Equation  (2.6)  implies  that  a  constant  spectral  distribution 
function  is  given  by  a  Dirac  delta  distribution  in  space,  and 
Eq.  (2.7)  implies  that  a  constant  function  in  the  spatial 
domain  has  a  spectral  distribution  function  given  by  the  Dirac 
delta  distribution.   This  interpretation  of  Eq.  (2.6)  and  Eq. 
(2.7)  follows  from  the  definition  of  Fourier  transform  pairs 
and  demonstrates  one  of  the  relationships  between  the  concept 
of  stochastic  processes  and  the  concept  of  Fourier  analysis. 
The  Fourier  transform  pair  for  a  real  spatial  function 
with  a  complex  spectral  density  function  will  be  denoted  by 

00 

?(±x)    =   /   F(^)    exp   i{±ic'x}dic  (2.8a) 

-oo 

F(±ic)    =   — i-^     /    C(x)    exp   i{+^-x}dx  (2.8b) 

(2U)^      -co 

In  place   of  the      Wiener-Khinchin   transform  pair    [5,26,88,131], 
the    Fourier  transform  pair   for  the    autocovariance    function 
and   the   spectral    density   function   for   a  homogeneous   process 
will  be    given  by   the    following: 

00 

Y(r)    =   /   S(t)    exp    i{it-r}dic  (2.9a) 

S(^)    =   — !— ^     /   Y(r)    exp   i{-t*r}dr  (2.9b) 

KK  (2tt)^    -co      cC 

The  norming  constant  required  for  the  inversion  of  a  Fourier 
transform  pair  will  always  be  associated  with  the  transform 
which  defines  the  spectral  density  function.   The  Dirac  delta 
distributions  defined  by  Eqs.  (2.6)  and  (2.7)  were  multiplied 
by  2Tr  in  order  that  their  Fourier  transform  would  be  propor- 
tional to  unity.   The  improper  integral  given  in  Eq.  (2.4)  is 
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often  associated  with  Bochner  (cf.  Wilks  [134])  or, if  the 
spectral  distribution  function  is  confined  to  an  interval 
{-7r,TT},  with  Herglotz  (cf.  Wilks  [134]). 

Comparing  the  stochastic  Fourier-Stieltj es  integral 
representation  of  a  random  function  given  by  Eq.  (2.1)  with 
the  Fourier  integral  representation  of  any  function  (i.e., 
both  deterministic  and  nondeterministic  functions)  which 
satisfies  the  requirements  for  the  existence  of  a  Fourier 
transform  pair  (cf.  Titchmarsh  [126]),  the  relationship 
between  the  random  measure,  dC(k), given  by  Eq.  (2.1), and  the 
Fourier  spectrum,  F(^), given  by  Eq.  (2.8a),  is  seen  to  be 

F(it)  =  dilll  (2.10) 

dk 

i.e.,    the    random  measure    of   the   process   must   have    a  density. 
Similarly,    comparing   the    spectral    representation   of   a  homo- 
geneous  process,    ds(k),    given  by   Khinchin's   Theorem  in  Eq. 
(2.4)    with   the    Fourier   transform   of  the    autocovariance    func- 
tion,   S(k),    given  by  Eq.    (2.9a),    the    relationship   again   im- 
plies   that   the    distribution   function  has    density 

S(^)    =   MM.  (2.11) 

?C  die 

The   major   distinction  between   a   random   function   represented 
by   a  Fourier-Stieltjes    integral    and   a  random  function   repre- 
sented by   a  Fourier   integral    is   the   existence   of  a   derivative 
of   the    distribution   function. 
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In   applications    to   random  surface   gravity  wave   problems, 
the    restriction   of  the    a   priori    assumption  that    the   process 
has    a  density   function   is    frequently   unimportant.      Kinsman 
[72;    Ch.    8.3]    discusses    the   mean   square    integrability   require- 
ment between   stochastic   and   Fourier   integrals    in   terms   of  the 
scales    of   the   process. 

The    shifting  property   of  the   Dirac   delta   distribution 
(cf.    Papoulis    [96])  plays    an   important   role    in   the    filtering 
problem  which    is    discussed   in   Chapter    3.      This    shifting  prop- 
erty  is    defined   in   the    space    domain  by 

(2TT)^6(x±r)    =    /  exp   i{it- (x±r)  }d^  (2.12) 

-  oo 

and   in    the  wave   number   domain  by 

&(t±tj    =  ^      f  exp   i{-x.  tf±it^)}dx  (2.13) 

°  (2tt)'^    -co  ° 

The    autocovariance    function   defined   above    is    one    of  the 
statistical   measures   which  will  be    required  to   compare    and   to 
characterize    a   random  process.      These    statistical   measures    are 
defined  by   an  ensemble    averaging   operator   denoted  by  the 
expectation   symbol,   e,    and   defined  by 

e{Z"}    =   /   C^p(Odc  (2.14) 

-oo 

where  p(C)  is  the  probability  density  of  the  random  variable 
Z.  The  variance,  a  ,  is  the  second  order  statistical  moment 
which  measures  the  dispersion  about  the  mean  and  is  computed 
from 
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ol   =  e{(Z-f{Z})2}  =  e{2^}-e'^{Z}  (2.15) 

The  standard  deviation  is  the  positive  square  root  of  the 
variance  as  defined  by 

STD  DEV  =  l^  (2.16) 

The  autocovariance  function  given  by  Eq.  (2.3)  for  a  random 
process  with  a  zero  mean  evaluated  for  zero  spatial  lag  r  is 
equivalent  to  the  variance. 

The  trivariance  function  is  a  third  order  statistical 
moment  defined  by  the  expectation  operator  as 

Y(x,r,,r2)  =  £{?(x)-c(x+r,)-^(x+r7)}  (2.17) 

The  skewness  is  a  measure  of  the  asymmetry  of  the  random  sea 
surface  about  the  mean  normalized  by  the  cube  of  the  standard 
deviation ,  i.e . , 

SK  =  ^»^-^">>'>  (2.18) 

The  skewness  measure  may  be  evaluated  from  the  trivariance 
function  for  a  random  process  with  zero  mean  for  zero  spatial 

lags  r^,r2- 

Although  higher  order  polyspectra  measures  may  be  defined 
(cf.  Brillinger  and  Rosenblatt  [29,30]),  only  the  kurtosis 
measure  will  be  evaluated  herein.   The  kurtosis  is  the  fourth 
statistical  moment  about  the  mean  normalized  by  the  square  of 
the  variance  and  is  referenced  to  an  equivalent  normal  random 
process  with  the  same  mean  and  variance  by  reducing  the 
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quotient  by  3;  i.e., 

K  =  ^^i1-^)^^^)^^    .  3  (2.19) 


3.   Tri variance  Function 

The  statistical  moment  which  gives  the  best  measure  of 
the  second  order  perturbation  contributions  to  the  random  non- 
linear sea  process  through  a  measure  of  the  asymmetry  of  the 
random  realization  about  the  mean  is  the  third  statistical 
moment  or  skewness.   Formally,  the  assumptions  of  stationarity 
and  homogeneity  results  in  making  the  functions  from  correla- 
tion theory  (i.e.,  covariance,  trivariance,  etc.)  dependent 
only  on  the  differences  in  the  temporal  and  spatial  parameters 
of  its  argument.   Although  the  assumptions  of  stationarity 
and  homogeneity  are  used,  the  explicit  dependence  of  these 
correlation  functions  on  the  temporal  and  spatial  parameters 
t  and  X  will  be  retained  in  order  that  the  relationship  between 
random  measures  and  the  spectral  distribution  function  given 
by  Eq.  (2.5)  may  be  employed.   This  definition  for  the  spec- 
tral distribution  function  involves  the  Dirac  delta  distri- 
bution as  a  result  of  the  assumptions  of  stationarity  and 
homogeneity.   The  Dirac  delta  distribution  greatly  facilitates 
the  derivation  of  certain  results  obtained  in  later  sections. 

The  trivariance  function  for  both  time  and  space  param- 
eters involves  two  lags  and  may  be  represented  by 

y(   T^,T2;   r^,r2)  =  E{ri(t,x)-Ti(t+T^;x+r^)-n(t+T2;x+r2)}  (3.1) 
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I£  the  random  sea  realization,  n(t,x),  is  assumed  to  be 
composed  of  a  linear  Gaussian  estimate  plus  a  nonlinear,  non- 
Gaussian  second  order  perturbation  component,  the  random 
process  may  be  expanded  in  a  power  series  given  by 

n(t,x)  =  ^n(t,x)  +  2n(t,x)  (3.2) 

where  the  order  of  the  perturbation  power  series  parameter 
has  been  implicitly  absorbed  into  the  random  sea  function  and 
the  order  of  the  perturbation  parameter  is  indicated  by  the 
preceding  numerical  prefix.   Substituting  Eq.  (3.2)  into  Eq. 
(3.1)  and  using  the  linearity  of  the  expectation  operator,  e, 
yields  the  following: 


+  £{j^ri(t;x)-^nCt+Tj^;x+r^)-2ri(t  +  T2;x+r2^^'^ 
+  E{^n(t;x)-  ^r]  (t  +  T^;x+r^)  •  ^n  (t+T2  ;x+r2)}+ 

+  £:{2n  (t;x)  -^n  Ct+T-|^  ;x+r^)  •  ^n  (t+T2  ;x+r2) }+ 

+  ^{0(5n)}  (3.3) 

where  ff{0(rn)}  means  of  order  five  in  the  perturbation  param- 
eter. 

Assioming  that  the  random  measure  of  the  linear  Gaussian 
estimate  is  uniformly  distributed  between  {-7t,7t},  the  Fourier- 
Stieltjes  integral  for  the  linear  estimate  may  be  replaced 
by  the  density  of  the  random  measure  to  obtain  the  following: 
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00 

^n(t,x)    =    ff¥[^]    exp   i{at-?-x}dadic 

■'•  -oo 

00 

=    //|F[ic]  l-exp   i{at-ic.x}-exp   i{-i,(t)  }dodt  (3.4) 

-  00 

where    ijj(ic)    is    an   independent    random  phase   angle   uniformly- 
distributed  between    {-tt.-it}.      Papoulis    [96,97]    demonstrates 
that   a  random   function  with   this    distribution   tends    rapidly 
toward   a  Gaussian   random   function.      The    first   expectation 
operator   in   Eq.     (3.3)    is,    therefore,    a   third  order   statisti- 
cal moment   of  a   Gaussian   random  variable   and   is    of   third 
order   in   the  perturbation  parameter.      Appendix  A  demonstrates 
that   differentiating  the   moment    generating   function    for   a 
Gaussian   random  variable    three    times   yields    a   zero  product 
for   this    term.      The   vanishing   of   this    term   for   a  Gaussian 
random   function  may   also  be    seen  by   substituting   Eq.     (3.4) 
into  the    first   term  of  Eq.    (3. 3). with   the    following   results: 

ffC^nCtjX)  •-|^n(t+T-,^;x+r-,^)-^n(t+T2;x+r2)      = 

fdo^    fdo^    /da3///|^F(ic^)t.|2F(ic2)|-|3F(^3)|- 

-co  -co  -00  -00 

•exp   i{  (a-j^+a2  +  a3)t- (^^+1^2+1^3)  •  x}  • 

•exp   i{a2T.j^-ic2^r-j^+a3T2-f  3'r2}  • 

•£{exp   i{-\l)^(t^')-\l)2(t^)-\l)^(t^)]dt^dt2dt^  (3.5) 

Since   the   random  phase    angles,    i|j(^),    are    independent,    the 
expectation   operator  may  be    decomposed   into   the   product   of 
three   expectations,    i.e.. 
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£{exp  i{-ip^(^^)}}-E{exp  i{ -i|;2  (ic2)}  •£{exp  i{-^^(t^)} 

(3.6) 
each  of  which  has  zero  expectation  and  is,  therefore,  equiva- 
lent to  the  probabilistic  results  derived  in  Appendix  A. 
Other  examples  of  uniformly  distributed  random  variables  are 
given  by  Borgman  [24]. 

These  results  demonstrate  an  important  method  of  measur- 
ing a  random  process  which  is  Gaussian  only  in  its  first 
order  estimate.   The  lowest  order  nonvanishing  statistical 
measure  of  the  trivariance  function  given  by  Eq.  (3.3)  for  a 
random  process  with  a  Gaussian  linear  estimate  is  of  fourth 
order  in  the  perturbation  parameter  and  is  due  solely  to 
existence  of  the  second  order  nonlinear  correction  to  the  linear 
Gaussian  estimate.   Moreover,  this  lowest  order  nonvanishing 
measure  is  closed  at  second  order  in  the  perturbation  param- 
eter as  no  other  higher  perturbation  orders  may  form  third 
order  statistical  moments  which  are  of  fourth  order  in  the 
perturbation  parameter.   In  order  to  interpret  further  the 
effects  of  the  second  order  correction  on  the  trivariance 
function,  an  explicit  form  of  this  correction  is  required  and 
will  be  computed  in  the  following  section. 
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4.      Nonlinear  Wave-Wave    Interactions    Correct 
to   Second  Perturbation   Order 

An  excellent    formulation   of  the   boundary  value   problem 
for  surface    gravity  waves    is    given  by   Wehausen    [127;    pp.    447- 
455]    and   this    formulation   is    equally  valid   for   random  surface 
gravity  waves.      As   noted   in   Section    1,    the    convergence    of 
random   functions,    derivatives    and  integrals   used   in   the    random 
boundary  value  problem  will  be    in  mean   square.      In   addition, 
Fourier-Stieltjes    integrals   which   require  no   a   priori    assump- 
tion  regarding   the   existence   of  a   density   function  will  be 
utilized. (cf.    Bharuchua-Reid    [13]). 

The    fluid   is    assumed   to  be    inviscid  and  incompressible. 
The    fluid  motion   is    assumed   to  be    irrotational    and   three- 
dimensional.      The    x,y   coordinate    axes    lie    in   a  horizontal 
plane    at   the   mean   sea   level,    z=0 ,   with   the    z   axis   positive 
upwards.      A   random  velocity  potential,    $,    exists    such   that 

u(x,y,z,t)    =    -^$  (4.1) 

For  notational  compactness,  the  horizontal  x,y  plane  will  be 
denoted  by 

X  =  xe^  +  yCy  (4.2) 

The  governing  partial  differential  field  equation  for  the 
random  velocity  potential,  $,  is  the  equation  of  continuity 
for  an  irrotational  fluid  flow; 

A$  =  v^$  =  0  ;   |x|«»,-h<z<n(x,t) ,t>0  (4.3) 
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where  the  three-dimensional  gradient  operator,  V(-),  is  given 
by 

'^'^-'-^K^V-^^'-^K  (4.4) 

and  e^    are    the    orthogonal   unit  vectors    in  the    three-dimensional 
cartesian    coordinate    system.      The    finite    depth   fluid   domain 
is    assumed  to  be   unbounded   in   the  horizontal   dimensions    so 
that   the  boundary   conditions    required   at   the   horizontal   boun- 
daries   located   at    infinity   are    those   which  preclude   partial 
standing  waves;    and  the   boundary   conditions    required  to  quan- 
tify  the    unknown   coefficients    or  eigenvalues    are   prescribed 
at   the    free   surface    and  horizontal  bottom  boundaries.      The 
fluid  is    assumed  to   occupy   a   finite    depth,    h,    in   the    lower 
half  plane    and   to   satisfy   the    following  no   flow  bottom  boun- 
dary  condition    (BBC)    across    the   horizontal    impermeable   bottom: 

BBC:      $^    =    0    ;      z    =    -h,|xi<   ~,t>0  (4.5) 

The  fluid  surface,  ri(t,x),  is  unconstrained  and  a  kinematic 
free  surface  boundary  condition  (KFSBC)  is  required  to  pre- 
scribe mathematically  the  continuity  of  fluid  motion  at  the 
interface  and  to  ensure  physically  that  no  fluid  particles 
are  convected  across  this  interface;  i.e., 

KFSBC:   n^-2^^-2^^''*z  "  °'   ^  "  n  (x,  t)  ,  |  x|  <»,  t>0        (4.6) 

where  the  two-dimensional  horizontal  gradient  operator, 
2^(0 »  is  given  by 


35 

An  additional  dynamic  free  surface  boundary  constraint  is 
required  to  ensure  that  all  stresses  acting  along  the  free 
surface  are  continuous.   In  the  absence  of  surface  tension 
and  for  an  atmospheric  pressure  assumed  to  be  zero,  the  dynamic 
free  surface  boundary  condition  (DFSBC)  may  be  obtained  via 
the  Bernoulli  equation,  i.e., 

-^^*gr\+j\'^^\^   =  Q(t);   z  =  n(x,t),|x|«»,t>0  (4.8) 

The  total  derivative  of  the  Bernoulli  equation  (cf.  Phillips 
[100],  p.  23)  combined  with  the  KFSBC  eliminates  the  unknown 
scalar  free  surface  elevation,  n(t,x),  from  the  DFSBC  boundary 
condition  to  yield  the  following  combined  free  surface  boun- 
dary condition  (CFSBC) : 

[jt  -   ^*-^](^t"S^'rl^'^I^^Q^^^)  =  °'   ^  =  n(x,t),|xl«»,t>_0 

(4.9) 
which  may  be  expanded  to  the  following  form: 

\t^g^2.^^t'\jt    "  7^  •  ^]l^^l^  =  0;   z  =  n(x,t),|x|<-,t>0 

(4.10) 
In  order  to  avoid  evaluating  the  CFSBC  at  the  unknown  free 
surface  elevation,  n(x,t),  and  in  order  to  utilize  a  separable 
coordinate  system  to  solve  the  governing  field  equation 
[Eq.  (4.3)],  the  CFSBC  is  expanded  in  a  Maclaurin  series 
about  the  still  water  level,  z=0,  i.e.. 
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v=o''h^  Kt^sVQt-t|t  -    1^    -^llv^l']  =  0;   z  =  o,|J|<»,t>o 

(4.11) 
In  order  to  obtain  approximate  solutions  to  Eq.  (4.3) 
subject  to  the  nonlinear  CFSBC,  the  nonlinearities  are  assumed 
to  impose  small  perturbations  on  an  initial  linear  approxima- 
tion.  The  measure  of  the  perturbation  is  assumed  to  be  pro- 
portional to  the  wave  slope  and  will  be  used  as  the  measure 
of  the  perturbation  ordering  parameter.   However,  it  is  not 
a  priori    clear  that  the  ordering  of  the  perturbations  by  a 
parameter  is  valid  for  either  a  continuous  or  discrete  spec- 
trum of  random  surface  gravity  waves.   Rellich  [108]  has 
devoted  considerable  research  into  providing  mathematical 
rigor  to  problems  of  perturbing  spectra  in  Hilbert  space  and 
has  established  a  sufficient  set  of  proofs  to  ensure  the 
validity  of  the  perturbation  concept.   Friedrichs  [50]  con- 
tinued the  work  initiated  by  Rellich  and  contributed  further 
to  the  completeness  of  the  proofs.   The  proofs  which  are 
required  to  establish  the  existence  of  a  perturbation  order- 
ing parameter  for  a  continuous  spectrum  in  Hilbert  space  are 
too  complex  to  be  covered  here.   The  material  given  by  Rellich 
[108]  and  Friedrichs  [50]  and  summarized  by  Dunford  and 
Schwartz  [43]  establish  the  validity  of  a  perturbation  order- 
ing parameter  for  a  continuous  spectrum  and  will  be  assumed 
without  proof  in  the  following  development. 

The  random  velocity  potential,  random  sea  surface  reali- 
zation, and  the  Bernoulli  constant  are  assumed  to  be  expand- 
able in  the  following  power  series  with  the  perturbation 
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ordering  parameter   implicitly   absorbed   into   the    functional 
notation: 

oo 

$(x,z,t)    =      I     ^$(x,z,t)  ^  (4.12a) 

v  =  l 

oo 

Ti(x,t)    =     I     ^n(x,t)  (4.12b) 

v=l 

Q(t)    =      E        Q(t)  (4.12c) 

v=l    ^ 

Substituting  Eqs .     (4.12a,b,c)    into   the   equations    of  the 
boundary   value   problem  and  equating  equal   orders   of  v,    the 
following   set   of   linear  boundary   value   problems    is    obtained: 

V  =    1:      V^$   =   0;      -h<z<0,|x|<~  (4.13a) 

^$^    =   0;      z    =    -h,  |x|<oo  (4.13b) 

l^tt'^Si^z'^l^t   "   °'      ^   "   O,|x|<=o,t>0  (4.13c) 

^n    =   I  ^i^t'^lQ^'    ^   "   0,|x|<-,t>0  (4.13d) 

V  =    2:      72$   =   0;    -h<z<0 ,  |  x|  <oo  (4.14a) 

2$^   =   0;    z    =    -h,  |x|<oo  (4.14b) 

2^t^g2^z^2Qt   =   It    IVI'-I^'Iy  ^l*tt^gl^>' 

z   =   0,  |x|  <oo,t>^0  (4.14c) 


1    .  I^^l'        1^ 

2^    =   1^2^-- T— "   i-l*tz^2Q>' 


z   =   0,  lx|<oo,t>0  (4.14d) 
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A  solution  which   satisfies    exactly  Eqs.     (4.13a,b)    is 
given  by   the    following   Fourier-Stielt jes    integral: 

^H^,z,t)    =    -i  7  g   .    coshClkKh.z)}    .  i(-^5^}.  A[lt,t] 

cosh{|lc|h}  1    ^    '    ■■ 

(4.15) 
where    d^A[k,t]    is    the    random  measure    of  the   velocity  poten- 
tial.     Equation    (4.15)    requires   no   a   priori    assumption   regard- 
ing either  the   existence    of   a   density   function   or  the   temporal 
dependence    of  the    random  velocity  potential. 

Absorbing   the    temporal    derivative    of  the    Bernoulli    con- 
stant  in   Eq.     (4.13c)    into   the    random  measure,    d,A[ic,t],    by 
Tequiring  that    the    random  process    maintain   a   zero   mean,    yields 
the    following   initial   value   problem   for  the    time    dependence    of 
the   random  measure: 

-i  /    g-exp   i{-k-x}-[^  +   g|k|tanh{|klh}]d,A[ic,t]    =   0      (4.16) 

The    deterministic   Green's    function  which   satisfies    this    ini- 
tial  value   problem  was    first   derived  by   Finkelstein    [45]. 
Friedman   and   Shinbrot    [47,48]    have   extended  his    results    to 
generalized  normal    domains    and  have    discussed  the    determinis- 
tic solutions    in  exhaustive    detail.      The   solutions    to   the 
random   initial   value   problem,    however,    are    less   perspicuous 
and   involve    an   irreversibility  paradox  with   regard   to   the 
statistics    of  the    initial    conditions.       In   order   to   avoid  the 
complications    attendant  with   these   initial   value    statistics, 
the    random  process  will   be    assumed   to  be  both   stationary   and 
ergodic.      These    assumptions    rer  '=r  the   solution  statistically 


39 

independent   of  any   initial   conditions.      A  discussion   of  the 
effect   of  these   assumptions   on  the   stochastic   integral   given 
by  Eq.    (4.15)    may  be    found   in   Kinsman    [72;    Ch.    8]. 

Equation    (4.16)    may  be    solved  by   the    application   of  the 
spectral   theory   of  random  operators    (cf.    Dunford   and   Schwartz 
[43]    or   Friedman    [49]).      The   integrand   in  Eq.     (4.16)    may  be 
written   in   spectral    operator  notation  by 

L[ic,t]-{djA}    =   0  (4.17) 

where   the    linear  spectral    operator,    l,    is    given  by 

i[k,t]    =   ^  +    g|k|tanh{|klh}  (4.18) 

at 

Equation  (4.17)  is  satisfied  either  by  the  trivial  solution 

d^A[ic,t]  E  0        V  it,t  (4.19) 

or  by   the    zeroes    of  the    spectral    operator 

L[t,t]    =0  (4.20) 

Equation    (4.17)    requires    that   the   spectrum  of  the   random 
measure,    d-,A[k,t],    vanish   identically   everywhere   except    along 
the  hypersurface   given  by   the    zeroes    of  the    spectral    opera- 
tor.     An   integral  which   satisfies   Eq.     (4.17)    may  be    repre- 
sented by   the    following: 

d^A[k,t]    =   /   exp   i{at}d^A[t,a]  (4.21) 

-00 

where  the  eigenvalues,  t,    are  determined  by  the  zeroes  of  the 
spectral  operator,  l;    viz. 
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a^  =  g|i^|tanh{|it|h}  (4.22) 

Equation  (4.22)  establishes  a  relationship  between  the  sto- 
chastic integral  given  by  Eq.  (4.15)  and  the  hydrodynamics 
of  the  boundary  value  problem.   The  hydrodynamic  boundary 
condition  requires  that  the  spectral  representation  of  the 
random  measure  vanish  everywhere  in  the  k,a  space  except  along 
the   curve    given  by  Eq.  (4.22).   This  restriction  on  the 
spectrum  of  d-,A[k,a]  may  be  given  by  the  following  Dirac  delta 
distribution  (cf.  Kinsman  [72],  Ch.  7.3). 

d^A[it(a)]  =  d^A[i^, a']6(a'-a)  (4.23) 

Alternate  methods  for  computing  the  eigenvalues  given 
by  Eq.  (4.22)  and  the  spectrum  of  the  random  measure  given  by 
Eq.  (4.2  3)  for  general  random  differential  equations  are 
given  by  Boyce  [26].   Methods  which  are  specifically  appli- 
cable to  random  waves  in  a  dispersive  media  are  given  by 
Frisch  [51].   For  the  random  surface  gravity  wave  problem 
formulated  above  the  fluid  media  is  assumed  to  be  nonrandom 
and  the  eigenvalues  computed  by  Eq.  (4.22)  are  deterministic. 
Additional  methods  for  applying  the  CFSBC  given  by  Eqs. 
(4.13c)  and  (4.14c)  are  given  by  Phillips  [100]  who  has 
pointed  out  the  similarity  of  the  CFSBC  to  the  equation  of 
an  undamped  harmonic  oscillator. 

Substituting  Eq.  (4.21)  into  the  general  solution  given 
by  Eq.  (4.15)  and  using  Eq.  (4.23)  gives  the  following  random 
velocity  potential: 
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,$(x,z,t)    =    -i  7g   •    cosh{|gUh.z)}    .    ^^p   i{at-it.5}H   ;i[ir(a)] 

^  -00  COSh{|ic|h}  1 

(4.24) 
where  the  reality  o£  the  random  velocity  potential  requires 
the  following  complex  conjugate  relationship: 

djA[-it(a)]  =  d^A*[i^Ca)]  (4.25) 

The   random  sea   surface    realization   is    determined   from  Eq. 
(4.13d),    i.e. , 

oo 

^n(x,t)    =   /  a-exp   i{at-ic-x}d^A[ic(a)]  (4.26) 

The  spectral  representation  of  the  random  measure  of  the 
random  velocity  potential  may  be  redefined  in  terms  of  the 
spectral  representation  of  the  random  sea  surface  elevation 
by  the  following: 

d'iF[t(a)]    =  a-djA[^(a)]  (4.27) 

In   terms    of   the    spectral    representation   of  the    random   sea 
surface   measure  ,  ^F[ic  (a)  ]  ,    the   random  velocity  potential    and 
random  sea   surface    realization  become    the   following: 

,Hx,z,t)    =    -i   7  f   .    cosh{|itKh.z)}    .    ^^p   i{at-it.^}d,F[it(a)] 
cosh{|klh}  ^ 

(4.28a) 

oo 

,n(x,t)    =   /  exp   i{at-ic-x}d,F[ic(a)]  (4.28b) 

-00 

The  functional  form  of  the  solution  to  the  second  set  of 
equations  which  are  correct  to  second  order  in  the  perturba- 
tion parameter  and  which  are  given  by  Eqs.  (4. 14a,b ,c,d)  is 
assumed  to  be  similar  to  the  first  order  solution  given  by 
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Eq.    (4.28a).      Denoting   the    frequency,    o),    and  the    separation 
constant,    kCw),    in    a   general   form  which  are  to  be   evaluated 
through   the   boundary   conditions,    a   solution  which   satisfies 
exactly  Eqs.     (4.14a,b)    is    given  by   the    following   random   inte- 
gral: 

2nx,z,t)    =    -i    7   cosh{|?Kh.z)}    .  i{a)t-?-5^}.d,s[?Ca))] 

^  -<»        cosh{|<|h}  2         ^    ^J 

(4.29) 

As   noted  previously,    Eq.     (4.14c)    is  equivalent   to   the 
equation    for   a   forced,    undamped  harmonic  oscillator.      This 
equation  has  been   studied  extensively   in  the   theory   of  spec- 
tral  operators    and  may  be  written   in   the  following   spectral 
operator   form: 

00 

-i    /  exp    i{a3t-K-x}-L[(jj,i<(co)]  •d2B[ic(a))  ]    = 

00 

=    //   exp   i{  (0^+02)  t-(.t^+t^)-x}-N[o^, 02',   ^1(0-^) '^2^'^  2^' 

'd^F[t^(o-^)]d2F[t2ia2)]  (4.30) 

where  L[a3,K(a))]  is  a  linear  spectral  operator  and 
iv[a,  ,a2  ;k^  (a.j^)k2  (02)  ]  is  the  spectral  operator  for  the  non- 
linear derivatives  given  by  the  right  side  of  Eq.  (4.14c). 
The  general  functional  form  of  the  separation  constant, 
K  (o))  ,  and  frequency,  u,  may  now  be  quantified  by  requiring 
equivalence  of  the  temporal  and  spatial  fionctions  in  Eq. 
(4.30).   This  equivalence  requires  a  change  of  variables  in 
the  integral  for  the  second  order  random  velocity  potential. 
Hildebrand  [64]  or  Sokolnikoff  and  Redheffer  [115]  give  the 
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requirements  for  the  application  of  the  following  change  of 
variables  equation: 


/  ax,y)dxdy  =  /  C{f(u,v),g(u,v)}|JC|^)|dudv 
R  R  "'^ 


(4.31) 


where  |J(')I  is  the  absolute  value  of  the  Jacobian  of  the 
change  of  variables  and  R,i?  are  the  regions  of  integration, 
Introduce  the  following  change  of  variables  in  Eq.  (4.30): 


<   =  ^i(o^)    +  ^2  (^^2^ 


which  have    the    following  nonvanishing   Jacobian; 


(4.32a) 
(4.32b) 


K  ,bi 


k^(a^) ,k^io^) 


3k 
9ki 

3k 
91^2 

8co 
9a, 

da. 

9a) 

da  2 

dic2 
,-th   , 

=  %-  %  (^-"^ 


where  C™  is  the  group  velocity  of  the  i    component.   The 
*^i 

CFSBC  now  becomes 


■i  //  exp  i{  (a^+a2)t- (icj^+ic2)  "x} 


^k^(a^),k2(a2)J 
•L[a^,a2;i^^,it2]-d2s[]^^(a^),ic2((J2^1  = 

00 

=  //  exp  i{  (a-j^+a2)t- (^^+ic2) -x}- 
'N[a^ya2;t^,t^]d^F[t^(a^)]d^F(.t^(a2)]  (4.34) 

which  may  be  solved  for  the  spectral  representation  of  the 
random  velocity  potential  correct  to  second  order  to  give  the 
following: 
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'    '   '    '   '     iJ^:^ '-^ >|-L[a,,a2;k,(a,),L(a,)] 

k^(a^),k2(a2)      ^   ^   1   1  '  2^  Z^J 

(4.35) 
provided  that  neither  the  Jacobian,  |J(*-)I,  nor  the  linear 
spectral  operator,  l[-],  are  singular.   The  diagonal  terms 
of  the  Jacobian  are  Stokian  nondispersive  waves  and  the 
Jacobian,  therefore,  may  be  seen  from  Eq.  (4.33)  to  be  non- 
singular. 

Phillips  [101]  and  Hasselmann  [54]  have  demonstrated  by 
separate  proofs  that  the  linear  dispersion  equation  given  by 
Eq.  (4.22)  is  convex  toward  the  k  axis  in  the  a ,k  plane  and, 
therefore,  no  zeroes  may  occur  in  l[-].   Dunford  and  Schwartz 
[43],  Friedman  [49],  and  Sveshnikov  [119]  give  polynomial 
expansions  for  the  ratio  of  spectral  operators  equivalent  to 
those  given  in  Eq.  (4.35)  whenever  zeroes  occur  in  the  spec- 
tral operator  l[-].   In  addition,  since  Eq.  (4.30)  is  an 
inhomogeneous  boundary  condition  resulting  from  the  pertur- 
bation expansions,  the  homogeneous  solutions  to  Eq.  (4.30) 
are  identical  to  those  given  by  Eq.  (4.17).   Therefore,  since 
the  spectral  operator  in  Eq.  (4.35)  is  nonsingular  and  since 
the  homogeneous  solutions  which  correspond  to  the  zeroes  in 
l[']    are  identical  to  the  linear  solution  and  limited  by  the 
Dirac  delta  distribution  given  in  Eq.  (4.23),  no  "free  wave" 
solutions  may  exist  at  second  order. 

Substituting  the  appropriate  derivatives  into  the  two 
spectral  operators  in  Eq.  (4.35)  yields  the  following  spec- 
tral response  function  for  the  random  second  order  velocity 
potential: 
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N[0^,0^S^,t^] 


L[0 


I'^Z'^l'^zl 


2Ca;^^a2^-[g^i^l-^2"^1^2^'''l^2'^''l^^^^g^^^l''2^^2^1^ 


2a^a2^(^'^l"'^2^    ^^l^l''^2l  '  tanh  [  |  kj^+k2  |h] } 

(4.36) 

Making  the  change  of  variables  given  by  Eqs .  (4.32a,b)  with 

the  Jacobian  given  by  Eq.  (4.33)  results  in  the  following 

equation  for  the  random  second  order  velocity  potential: 


(x,z,t) 


<»  cosh{|k^+k2l  (h+z)} 

:   -i   //   _ _ 

cosh{ |k^+k2 |h} 
■exp   i{ (0^+02) t- (^^+^2^ *^^* 


^lf^l^'^2'^^2^ 


ci2s[k-^(a^),k2(a2)] 


(4.37) 


Substitution  for  the  random  spectral  representation, 
d2B[ic,  (a, )  ,^2  (^^2-^  ■' '  Siven  by  Eq.  (4.35)  and  for  the  spectral 
response  function  given  by  Eq.  (4.36)  yields  the  following 
expression  for  the  random  second  order  velocity  potential: 


2$(x,z,t) 


i  // 


cosh{  11^^+1^2  1  (h  +  z)} 


cosh{  |k-|^+k2  |h} 
•exp  i{  io  ^+a  2)  t-  (tj^+t^)  'x} '  d[o  ^,0  2;t^(a  -^)  ,^2^0  2)  ]  ' 

•d^F[5^(a^)]-d2F[it2(a2)]  (4.38) 

where  the  nonlinear  interaction  kernel,  ^[a,  ,02  ;k,  (a^)  ,k2  ((^2)  ]  , 
is  equivalent  to  the  spectral  response  function  given  by  Eq. 
(4.36). 
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The  second  order  correction  for  the  random  sea  surface 
realization  may  now  be  evaluated  by  substituting  into  Eq. 
(4.14d)  the  proper  derivatives  of  Eqs.  (4.28a,b)  and  Eq. 
(4.38).   Evaluating,  as  before,  the  second  order  Bernoulli 
constant,  -Q*  in  order  to  maintain  a  zero  mean  second  order 
random  process  yields  the  following  expression  for  the  second 
order  random  sea  realization: 

2n(x,t)  =   2^  ^J"   exp  i{(o^  +  a^)t-(^^+^^)'x}-H[o^,a^'X^,t^]- 

where    the   nonlinear   interaction  kernel,    ^[a,  ,a-;ic,  jic^]  ,    is 
the    spectral    response    function   determined  by   the    derivatives 
of   the    right   side   of  Eq.     (4.14d)    and   is    given  by 

•Z)[a^, 02  ;i^^(a^)  ,^2(^02)]-      g  g   ^   •"^i^z^'^l'^^Z  (4.40) 

Note  that  the  second  order  random  sea  realization  given  by 
Eq.  (4.39)  is  inversely  proportional  to  the  gravitational 
constant,  g. 

The  nonlinear  interaction  kernel  given  by  Eq.  (4.40)  is 
expressed  in  a  symmetric  form.   Different  forms  may  be  ob- 
tained when  the  appropriate  derivatives  are  substituted  into 
the  inhomogeneous  forcing  term  on  the  right  side  of  Eq. 
(4.14d).   Wiener  [132]  has  pointed  out  that  any  kernel  oper- 
ating on  the  product  of  identical  random  functions  having 
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symmetric  domains  o£  integrations  may  be  rendered  symmetric 
by  reversing  the  arguments  in  the  kernel,  adding  the  two 
kernels  and  dividing  by  two.   This  procedure  has  been  used 
to  obtain  the  symmetric  kernels  used  in  Eqs .  (4.36)  and  (4.40). 

The  equation  for  the  nonlinear  sea  surface  realization 
correct  to  second  order  is  obtained  by  adding  Eqs.  (4.28b) 
and  (4.39) ,  i.e. , 

oo 

n(x,t)  =  ^Ti(x,t)+2n  (x,t)  =  /  exp  i{at-ic-x}-dF[^(a)]  + 

+  J-   //■ff{a^,a2;ic^(a^)  ,1^2(^2^^' 

•exp  iiia  ^+0  ^)t-  (t^+t^)  -x]  •  d^F[^^(a  ^)  ]  •  d^F  ["^^ia  ^)  ] 

(4.41) 
Upon  making  the  following  change  of  variables  in  the  integral 
for  the  second  order  random  sea, 

(4.42) 

the   expression   for  the    random  sea   realization  becomes 

00 

n(x,t)    -   /    (exp   i{-ic(a)-x}-dF[i^(a)]    + 
-  00 

^   Jg     ^   «{c^l,c^-(Ji;i^l((7i),it2(a-a^)}- 

•exp   i{  -  [i^-^  ((7^)  +i^2  (^^"^1^  ]  '^^  '  ^i^Li^i  i^i)  ]  •  d2F  [^2  (^'^1^  1  ^  * 

•exp   i{at}  (4.43) 

Comparing  Eq.     (4.43)   with  Eq.    (2.1),    the    random  non- 
linear  sea  surface    correct   to   second  order   is    seen   to  be 
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oo 

n(x,t)  =  /  exp  i{at}-dx[^io) ,x]  (4.44) 

-  00 

where 

dx[^(a)  ,x]  =  exp  i{-ic-x}  •  dF[ic  (a)  ] 

1   "  ^      ^ 

•exp  i{-[it^+i^]-x}-d^F[it^(a-L)]-d2F[lt2(a-a^)]  (4.45) 

The  important  point  to  note  about  this  final  expression 
for  the  random  sea  surface  realization  is  that  for  a  fixed 
set  of  horizontal  spatial  coordinates,  the  quadratic  contri- 
butions occur  at  the  same  frequencies  as  those  of  the  linear 
Gaussian  realization.   This  is  a  key  point  in  the  application 
of  the  fast  Fourier  transform  and  the  resulting  simplification 
of  the  computation  of  nonlinear  random  realizations.   Note 
also  that  the  second  order  contributions  to  the  random  mea- 
sure, dx[k(a),x],  which  are  given  by  the  integral  term  in 
Eq.  (4.44),  are  a  function  of  the  sums  and  differences  in 
wave  numbers,  k(a).   Consequently,  all  of  the  second  order 
contributions  which  result  from  this  convolution  integral  are 
not  phase  locked  or  nondispersive  with  respect  to  the  linear 
Gaussian  estimate  at  that  same  frequency.   The  simple  example 
for  deep  water  waves  given  in  Section  2  of  Chapter  1  illus- 
trated the  manner  in  which  these  second  order  contributions 
are  nondispersive  with  respect  to  the  linear  Gaussian  spectrum 
and  the  resulting  linear  profile. 
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5.      The   Bispectrum 

The  bispectrum  is    the   Fourier   transform  of  the    trivari- 
ance    function   and   represents    the   spectral   density   of  the    con- 
tributions   to   the   mean   cube   of   the    random  sea  realization, 
n(x,t),by   a  product   of   three   Fourier   densities   whose    resultant 
frequency   equals    zero    (cf.    Hasselmann   et^  al.    [57]).      This 
Fourier   transform  pair  was    shown   in   Section   3   to   yield   a  mea- 
sure  of  the    second  order   contributions    to   the    random  sea   reali- 
zation which   is    closed   to   the    second  order   in  perturbation 
parameter.      Since   much   useful    information   regarding   stationary, 
homogeneous    random  processes   may  be    obtained   from  this    trans- 
form pair,    a  brief  discussion   of  the  bispectrum   is    developed 
in   this    section. 

The    lowest   order  nonvanishing   third   order   statistical 
moment    for   a   linear  Gaussian   estimate  was    computed   in   Eq. 
(3.3)    and  may  now  be    decomposed  explicitly   into   first    and 
second   order  perturbation   contributions    using   the    results    of 
Section   4: 

yCt,  ,T^;r,  ,r^)    = 

+H[o2,o^;t2,t^]'exp    i{  (02  +  03)1^+04X2  -  (^2"'^3^^1 '^4  "^2  ^"^ 
.+ff[aj^,a2;ic-|^,ic2] -exp    ^^(^;^'^i-^o ^T^-t^-T^-t^-T^} 


ffff 


•exp    iiio^+a  2*0  ^+a  ^)t-  i'^^+t^+^^+t^)  -1}  ' 
•E{d^F[^;L^a^)]-d2F[it2(a2)]-d3F[^3(a3)]-d4F[^4(a4)]} 
*0{^T^}  (5.1) 
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This  is  a  fourth  order  moment  of  the  Gaussian  random 
measure,  dF[k(a)],  which  may  be  factored  as  shown  in  Appendix 
A  into  the  sum  of  three  product  pairs  of  all  possible  permu- 
tations of  the  four  Gaussian  variates.   If  the  process  is 
assumed  to  be  spatially  homogeneous,  the  spectral  distribution 
function  may  be  expressed  by  Eq.  (2.5)  in  which  the  arguments 
of  the  Dirac  delta  distributions  are  all  possible  permutations 
of  wave  numbers  which  result  from  the  products  of  the  decom- 
position of  a  fourth  order  Gaussian  random  measure.   The  fre- 
quency dependence  of  the  wave  number  given  by  the  linear 
dispersion  equation  in  Eq.  (4.22)  also  results  in  a  temporal- 
ly stationary  process  when  the  integrations  over  the  Dirac 
delta  distributions  are  performed.   This  fact  has  already 
been  used  to  eliminate  the  frequency  dependence  from  the  ran- 
dom spectral  measure  in  Section  4.   It  can  be  shown,  finally, 
that  the  term  of  lowest  order  in  the  perturbation  parameter 
which  is  nonvanishing  in  the  trivariance  function  for  homo- 
geneous, stationary  random  surface  gravity  waves  reduces  to: 

oo 

nnn  -<»  nn       nn 

•F[a]^,a2,T-j^,T2;^l,^2'^l'^2^^^1^^2       (5.2) 

where    the   nonlinear   interaction  kernel   F[ai  ,a2»Ti  ,T2 '»^-i  »^2 » 
r,  jT-]    includes   nine   permutations   of  the   previously   derived 
nonlinear   interaction  kernel,    // [a-,  ,a2  ;k,  ,1^2  ] ,    and   the    complex 
unit   vectors   of  the    time    and  space    lags,    t    and   r. 
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Formally,  the  double  Fourier  transform  of  the  trivari- 
ance  function  is  the  bispectrum  and  may  be  expressed  by 

exp  iC-o),  T, -W2T2+k,  •  r, +k2- r2}dr,  dr2dT-|^dT2         (5.3) 

Hasselmann  et  al.  [5  7]  have  given  the  symmetry  proper- 
ties of  the  Fourier  transform  pair  for  the  trivariance  func- 
tion and  bispectrum  and  have  investigated  the  bispectrum  of 
pressure  records  from  surface  gravity  waves  to  determine  the 
directional  spread  of  the  wave  spectra.   Brillinger  and 
Rosenblatt  [29,30]  give  smoothing  functions  and  other  useful 
relations  for  these  general  higher  order  polyspectra. 
MacDonald  [82],  Hinch  and  Clay  [66],  and  Haubrich  [63] 
derive  useful  expressions  for  determining  bispectral  esti- 
mates via  the  Fast  Fourier  Transform  algorithm.   Although 
much  useful  information  about  the  random  sea  process  appears 
to  be  discernible  from  bispectral  estimates,  only  the  coarse 
measure  of  the  trivariance  function  for  time  and  space  lags 
which  are  identically  zero  (viz.,    the  skewness)  will  be 
employed  in  this  study. 

Attempts  to  obtain  a  convolution  relationship  from  Eq. 
(5.2)  for  the  bispectral  estimates  given  in  Eq.  (5.3)  which 
would  be  similar  to  those  given  by  Tick  [123]  for  the  second 
order  spectra  did  not  prove  successful. 
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6 .   Second  Perturbation  Order  Spectra 

Although  the  autocovariance  function  for  random  sea 
surface  realizations  is  not  closed  at  the  second  perturba- 
tion order,  it  is  of  interest  to  compute  the  spectral  contri- 
bution which  results  from  the  second  perturbation  order  cor- 
rection to  the  random  linear  sea  realization. 

The  autocovariance  function  for  homogeneous,  stationary 
random  sea  correct  to  second  perturbation  order  is  given  by 

Y(T,r)  =  £:{[^n(t,x)+2ri(t,x)]- [^n(t+T;x+r)+2nCt+T;x+r)]} 

=  £{^Ti(t;x)^nCt+T;x+r)}  + 

+  £:{^nCt;x)  '211  (t+x  ;x+r)+2Ti  (t  ,x)  •  ^n  (t+x  ;x+r) }  + 

+  £:{2n  (t  ;x) '211  (t+T;x+r)} 

=  ////  exp  i{  (0-1^+02) t- (5^+ic2)  •x}-exp  1(02^-1^2'^^' 
-  00 

00 

+    ffffffff  exp    i{(a^+02+OT^+a^)t-(t^+t2+'^^+t^)'x}- 
-00 

•    exp   i{  (a3+a4)T-(lc3+54)-r}- 

'    E{d^F[t^,a^]'d2F[t2,a2]'d^F[t^,a^]-d^F[t^,a^]}    (6.1) 

since  the  expectation  of  the  product  of  the  first  and  second 
order  contributions  is  a  third  order  Gaussian  moment  of  inde- 
pendent random  variables  and  identically  equal  to  zero. 
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Using  the  identity  between  a  random  measure  and  the  spec- 
tral distribution  given  by  Eq.  (2.5)  and  the  decomposition  of 
a  Gaussian  fourth  order  moment  given  in  Appendix  A,  the 
following  spectral  representation  of  the  autocovariance  func- 
tion for  a  random  sea  correct  to  second  order  is  obtained. 

oo 

Y(T,r)    =   ////   exp    i{(a^+a^)t-(t^+t^)'x}'exTp   i{a2T-ic2 -x}  • 
nn 

•  S[i^2''^2^*'^'^^2''*^l^'^'^^2''^l^''^^2'^^l^'^2^^1    * 

oo 

+   ffffffff   exp    i{(a^  +  a2  +  a2+a4)t-(it^+ic2+k3+ic4)-x}- 

-  00 

•  exp    i{Ca3+a^)T-(ic3+ic^).r}- 

•  {//[a^,a2;k-,^Ca-j^),k2C^2^^'^^^3'^4'^3^'^3^'^4^^4^^' 

.    {S[it2'^2l''^^^2'"^l^"^^^2'"^l^*^t^4'^4l"'^^^4''^3^''^'^'^4''^3^ 

+   s  [1^4,0^] -5(1^4+1^^).  6  (a4+a^)-s  [1^3,03] -6  (1^3+1^2^  • 

•  6(03+02)}}    'i^i  •  ^1^2  •  dk3- dk.  •  d0,  •  d02- dQj- d0^  (6.2) 

Integrating   and  employing   the    shifting  property  of  the   Dirac 
delta   distribution,    the    following   expression   for   the    auto- 
covariance   is    obtained: 

00 
Y(T,r)    =    //  exp   i{-0,  T+k  •  r} 's  [-k-|  , -0-,  ]dk,  d0,  + 
nn 

00 

+    ////  exp   i{-(0-j^+02)T+(^-|^+k2)-r}-s[-^^,-0-|^]  •s[-ic2,-0  2] 

-00 

(continued) 
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•/f[-a^,-a2  ; -k^,-k2]  -exp   H- (0^  +  02)1+ (^^+^2)  •r}  + 

+if  [0-^,02  ;i^^,ic2]-H[-a2,-a^; -^2 '-^ll* 
•exp   i{-io2  +  o^')x+(t.2+t^')  -r} 


dk,dk2da^da2 


(6.3) 


The    synunetric  properties    of  both    the   spectral   distribution 
function,    5[k,a],    and   the   nonlinear   interaction  kernel, 
ff [0^,02  ;k-|^,k2]  ,    reduces   Eq.     (6.3)    to   the    following: 

00  00 

Y(T,r)    =    //   exp    i{-aT+^-r}-5[i^,a]d^da+    ////   s[ic,,a,]' 

rin  -°°  -00  -•■      -•■ 

•S[k2,cf2]  ■^['^p-c']L'^r'^l^  •^[a2,-a2;ic2,-^2^'^^1^^2^^1^'^2'' 

+  2////  exp   i{  -  (0^+02 )  T+  (ic^+i^2 )  •  r>  •  5  [i<^ , a^]  • 


.„2 


•5[k2,cr2]*^    [a^,a2;k,,k2]dk,dk2da,da. 


(6.4) 


Formally,    the   Fourier   transform  of  the    autocovariance    function 

yields    the    spectral    density   function.      Thus,    multiplying  both 

sides   of  Eq.     (6.4)    by   exp   i{-coT  +  K'r}    and   integrating   over 

time    and   space   yields 

00 
//   Y(T,r)    exp    i{ -a)T+i<- r}dTdr  = 
-<»     nn 

CO 

=    ////   5[^,a]'exp    i{- (a+a))T+(^+ic)  •r}d^dadTdr+ 
-00 
00 

+   ffffff  s[ic-j^,a^]  •s[^2'^2-' '^^"^l^'^l'^l' '^l-l  * 

•    H[a2  , -02  ;i^2, -ic2] -exp    i{ -cox+ic  •  r}dic,  d^2^'^i 'i'^2^^'^^'^ 

(continued) 
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°°  ->  ->-  2  ->■      -> 

-00 

•  exp   i{- (aj^+a2+w)T+(ic-j^+ic2  +  K)  •r}dk^dk2da^da2dTdr      (6.5) 

Inverting  the    order  of   integration   over   the  wave   number   and 
frequency   differentials  with   the    time    and   space    differentials 
results    in   integrals   equivalent   to   the    Dirac   delta   distribu- 
tion  given  by  Eqs.     (2.13)    for  the    spatial   integration   and   a 
similar   indentity    for   the    temporal    integration.      Substituting 
these    indentities    for   the    Dirac   delta   distribution,    the 
following   expression    is    obtained: 

//  Y(T,r)    exp   i{-uT+K' rldrdr   = 
-00     nn 

CO 

=    (2tt)^    //   s[i^,a]-6[i^+K)-5(a  +  u)dicda+ 

-co 

00 

+    (217)^   ////   s[k^,a^]-s[i^2'^2^'^^'^l''^l'^l'"^l^* 

•  H[a2,-a2  ;^2»  "^2-' '  '^  (^'^)'^  {y:')^^^^^Q ^^o ^■■<- 

+   2(27T)^   ////   s[it^,a^]-s[it2'^2^'^^t''l'''2'^l'^2^* 

•  6  (a^+a2+w) -6  (^^+ic2  +  <)dic^dk2da^da2  (6.6) 

Dividing  both   sides   by    (2it)      and   comparing   the   expression   for 
the    generalized   Fourier   transform  pair   for   the    spectral    den- 
sity  functions   yields    the    following: 

S[a),i<]    =   J  //  Y(T,r)*exp    i{ -oox+ic- r}dTdr 

nn  (2Tr)  nn 

The    first    integral   in   Eq.    (6.6)    gives 

,S[w,i<]    =   s[a),K]  (6.7) 

nn 
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which  is  the  spectral  density  for  the  linear  Gaussian  esti- 
mate.  The  second  integral  vanishes  identically  because  the 
argument  o£  the  Dirac  delta  distribution  is  not  in  the  domain 
o£  the  integral.   Integrating  the  third*  integral  and  using 
the  shifting  property  o£  the  Dirac  delta  distribution  gives 

00 

2S[k,w]  =  2//  s[k,  ,a,  ]  •s[K-ic,  ,a)-a,  ]  • 

•  H    [a,  ,w-a, ;k, ,<-k, ] • dk, da,  (6.8) 

Equation  (6.8)  is  the  important  result  first  given  by  Tick 
[123]  for  a  continuous  spectrum  of  deep  water  waves  and  later 
by  Longuet-Higgins  [79]  for  a  discrete  spectrum  of  deep  water 
waves  that  the  second  order  spectrum  which  results  only  from 
first  and  second  perturbation  order  contributions  is  a  con- 
volution of  products  of  all  first  order  spectral  densities 
whose  sums  and  differences  in  frequencies  contribute  to  the 
same  frequency  of  the  second  order  spectral  density.   Equation 
(6.8)  is  a  general  expression  for  the  second  order  spectral 
density  for  both  time  and  space  lags.   The  hydrodynamic  rela- 
tionship which  restricts  the  spectrum  of  the  random  measure 
to  those  values  in  the  a,k  plane  given  by  the  linear  disper- 
sion equation  has  not  been  invoked. 

The  simulation  of  a  time  series  for  a  random  sea  requires 
that  the  spectral  density  be  given  in  the  frequency  domain 
vice  the  wave  number  domain.   The  spectral  density  has  been 
given  in  the  wave  number  domain  as  a  result  of  solving  the 
spatial  boundary  value  problem  in  Section  4.   The  transforma- 
tion from  wave  number  space  to  frequency  space  will  be  con- 
sidered next. 


CHAPTER  3 
APPLICATIONS  OF  NONLINEAR  RANDOM  SEA  SIMULATIONS 

Section  1  transforms  the  directional  random  measure 
expressed  in  wave  number  space  for  the  linear  Gaussian  random 
sea  into  the  spectral  density  expressed  in  frequency  space. 
Section  2  discusses  the  fast  Fourier  transform  algorithm  and 
its  application  to  synthesizing  random  sea  realizations. 
Section  3  introduces  the  Bretschneider  spectrum,  evaluates  the 
two  constant  parameters  required  to  represent  measured  hurri- 
cane generated  wave  spectra,  and  discusses  the  application  of 
the  Phillips  equilibrium  spectrum  in  the  simulation  of  non- 
linear random  seas.   Section  4  briefly  describes  the  measured 
hurricane  generated  wave  and  pressure  force  realizations 
recorded  by  Wave  Project  II  during  Hurricane  Carla  in  the 
Gulf  of  Mexico.   Section  5  compares  the  distributions  and 
spectra  of  measured  sea  surface  realizations  from  four  hurri- 
cane generated  records  with  the  simulated  linear  and  nonlinear 
random  sea  surface  realizations.   Section  6  derives  the  equa- 
tions required  in  the  application  of  the  digital  linear  filter 
technique  to  predict  wave  pressure  forces  from  random  sea 
surface  realizations.   Section  7  compares  the  distributions 
and  spectra  of  the  measured  and  predicted  wave  pressure  forces 
at  the  55.3  feet  dynamometer  elevation  on  a  vertical  piling 
located  in  99  feet  of  water. 
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1.   Fourier  Series  Approximations 
of  the  Random  Measure 


Although  the  nonlinear  interaction  kernels  for  water 
of  finite  depth  computed  by  Tick  [125]  and  by  Hasselmann  [54] 
have  been  available  for  several  years,  no  successful  attempts 
to  employ  these  kernels  to  construct  second  order  nonlinear 
sea  simulations  have  been  published.   Tick  [125]  attributed 
this  to  the  computational  difficulties  involved  in  applying 
the  theory. 

The  availability  of  the  Fast  Fourier  Transform  (FFT) 
algorithm  [37]  now  makes  the  simulation  of  nonlinear  random 
seas  correct  to  second  order  in  the  perturbation  parameter 
relatively  simple.   Moreover,  the  retention  of  complex  phase 
angles  by  the  FFT  algorithm  makes  it  possible  to  compute  the 
spectral  densities  correct  to  second  order  entirely  in  the 
frequency  domain  and  to  make  one  inversion  to  the  time  domain 
for  the  final  sea  surface  realization.   Although  no  simultane- 
ous temporal  simulations  of  the  sea  surface  at  varying  hori- 
zontal spatial  locations  were  made  in  this  study,  the  exten- 
sion of  the  technique  developed  in  this  chapter  is  trivial  as 
will  become  obvious  later.   These  simultaneous  simulations  at 
varying  horizontal  spatial  locations  are  required  for  the 
dynamic  analysis  of  multilegged  pile-supported  platforms. 

Simulation  by  digital  computer  of  nonlinear  random  seas 
correct  to  second  order  using  Eq.  (4.43),  Chapter  2,  Section  4 
requires  that  the  exact  Fourier-Stieltjes  integrals  be  approxi- 
mated by  Fourier  series.   The  process  of  transforming  the 
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Fourier-Stieltjes    integrals    requires    first   that   the   spectral 
distribution    function  have    density: 


dF[k(a)]    =   F[k(a)]-dk 


(1.1) 


The  spectrum  of  the  random  measure,  F[k(a)],  is  in  wave  number 
space  as  a  result  of  solving  a  spatial  boundary  value  problem 
in  Chapter  2  and  a  transformation  to  frequency  space  is 
required  in  order  to  simulate  a  time  series.   The  transforma- 
tion from  rectangular  cartesian  coordinates  to  polar  coordi- 
nates is  given  by 


|k|  cos 


ky  =  |k|  sin 


with  the  following  change  of  variables  Jacobian; 


(1.2a) 
(1.2b) 


3(k,,V 


t3(|k|,6) 


cos    e  sin   e 

-|ic|    sin    e         |]c|    cos 


m 


(1.3) 


Substituting  this  change  of  variables  yields  the  following 
discrete  polar  spectral  representation: 


F[k(a)].dk  =  F[k,e]-k-dk.d6 


(1.4) 


Finally,  the  transformation  from  the  scalar  wave  number  space 
k,  to  wave  frequency  space,  a,  is  made  via  the  linear  dis- 
persion equation  to  obtain  the  following: 
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F[k,6]-k-dk.de  =  F[a,e]-  |J{|^}|-da-de  (1.5) 

where  the  Jacobian  of  the  change  of  variables  is  the  follow- 
ing: 

|J{|^}|    =    {C^}"^   =    {C    [1    +   2_[il^ ]}-!  (1.6) 

^"^  S  ^  sinh{2|it|h} 

The    determination   of   the    directional    dependency   of   a 
general   three-dimensional    random  process    is   not   trivial   even 
for   a   linear  Gaussian   estimate.      A  good   description   of  the 
theoretical    derivation    is    given  by   Longuet-Higgins    [78]    which 
is   based  on   the    theories    of  Raleigh    [105]    and  Rice    [109]. 
Additional    insight   into   the    directional   spectrum  problem  with 
some    specific   applications    is    given  by   Kinsman    [72;    Ch.    8]. 
Applications    to   directional  wave    arrays    are   given  by   Borgman 
[19]    and  Snyder   and  Smith    [121].      Snyder   and  Smith    [121]    out- 
line   an   extension   of  the    linear  directional    spectrum  to   a 
nonlinear   directional    spectrum   for   the    case    of   an   array   of 
wave    recorders.      Hasselmann   ejt    al_.     [57]    give    a  nonlinear   appli- 
cation   for   the    determination   of  the    directional    spread   from  a 
single   pressure    realization.      These    references    reflect    the 
complexities    involved   in   describing   the    directional    dependency 
of  the   spectral    representation   of  a   random  sea. 

For  the  simulation  of  a  time  series  of  a  random  sea,  the 
horizontal  spatial  coordinates  may  be  considered  as  fixed  and 
the  directional  dependency  may  be  expressed  implicitly  in  the 
random  measure  by   the    following   integral: 
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F[a]  =  /  k^I^-^  -exp  i{-k(cos  (e)x^  +  sin(e)y^)  }de      (1.7) 
o      g 

where  X, ,y,  are  fixed  values  of  the  horizontal  spatial  coordi- 
nates. 

The  approximation  of  the  stochastic  integral  representa- 
tion of  a  random  sea  realization  by  a  Fourier  series  requires 
that  the  discrete  stochastic  amplitude  of  the  approximating 
series  represent  the  equivalent  energy  content  as  that  con- 
tained in  the  exact  integral  representation  over  a  differen- 
tial interval,  da,  i.e., 

|F^Cn)|^  =  iF[nAa]l^-Aa  (1.8) 

Including  the  directional  dependency  implicitly  through  the 
complex  phase  of  the  stochastic  amplitude,  F,  (n)  ,  the  sto- 
chastic integrals  given  by  Eq.  (4.41)  in  Chapter  2,  Section  4, 
may  be  approximated  by  the  following  series: 


n  (t)  =   Z   F^(n)expi{nat}  + 
n=-oo 

00       oo 

+  1    E     E    F  (n)F  (m)H[m,n;k  ,k  ]• 
&  n=-«>  m=-<» 

•    exp   i{(n+m)at}  (1.9) 

The    simulation  by   digital    computer  of  the   random  process 
represented  by   the   stochastic   series    in  Eq.     (1.9)    requires    the 
evaluation   of  a   suitable   expression    for   the   stochastic   random 
function,      F^(n),    and   also   requires    that    this    function   repre- 
sent  a  Gaussian   random  process.      Descriptions   of  various 
stochastic   amplitudes   which  may  be   used  to   simulate    the    random 
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linear   sea   are    given  by   Birkhoff  and  Kotik    [14],    St.    Denis 
[117,118],    Le   Mehaute    [74],    Kinsman    [72],    inter   alios. 
Techniques    for   simulating  by  digital    computer   a   random 
Gaussian   realization   utilizing   these   stochastic   amplitudes 
or  deterministic   amplitudes    are    given  by  Mihram    [89]    and 
Borgman    [20,2  3]. 

The   technique    selected   for   simulating   linear  Gaussian 
random   seas    for   the   purpose   of   this    study   is    the    following: 

(1)  Generate    a   sequence   of  random  numbers   which 
are   uniformly   distributed  between    {0,1}. 

(2)  Multiply   the   uniformly   distributed   random  number 
sequence   by   2tt   in   order   to   obtain   a   random  phase 
angle   which   is   uniformly   distributed  between 
{0,2tt}. 

(3)  Form  a   complex  number  of  unit    amplitude   via 
Euler's    equation   using   the    uniformly   distributed 
random  phase    angle    as    the    argument. 

(4)  Multiply   the    complex  number   computed   in    (3)    by 
the    square    root   of  the    differential   energy   con- 
tent   given  by   the   Bretschneider  spectrum. 

Following   this    simulation  procedure,    the   stochastic   complex 
amplitude    in   Eq.     (1.9)    may  be   expressed  by   the    following   two- 
sided  Fourier   spectrum: 

F,(n)    =      ijs(n)    .    V-  •     [exp   i{-4;(n)}]  (1.10) 

^  I  BB  ^R 
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where  S(n)  is  the  value  o£  the  one-sided  Bretschneider  spec- 

BB 
trum  at  frequency  n-Zu/Tj^,  Tj^  is  the  length  of  simulated 

record,  and  ^(n)    is  the  random  phase  angle  generated  in  pro- 
cedure C2)  above. 

This  procedure  may  be  shown  to  be  equivalent  to  filter- 
ing Gaussian  white  noise.   Details  may  be  found  in  Bendat  [5], 
Davenport  and  Root  [38],  Laning  and  Battin  [73],  inter   alios. 

By  referencing  the  horizontal  spatial  coordinate  to  a 
fixed  location,  (x=0,  say),  the  complex  stochastic  amplitude 
for  the  random  sea  simulation  may  be  computed  by  initializing 
the  complex  coefficients  of  the  fast  Fourier  transform  with 
the  linear  stochastic  values  given  in  Eq.  (1.10)  plus  the 
second  order  contributions  given  by  the  double  summation  in 
Eq.  (1.9).   In  order  to  understand  this  method  better,  a  brief 
description  of  the  fast  Fourier  transform  algorithm  is  needed 
and  will  be  presented  next. 

2.   Fast  Fourier  Transform  (FFT) 

The  fast  Fourier  transform  algorithm  used  to  simulate 
the  random  nonlinear  sea  surface  profile  is  the  Subroutine 
NL0GN  written  by  Robinson  [110].   Since  most  F0RTRAN  compilers 
do  not  admit  zero  or  negative  subscripts  for  arrays,  the  nota- 
tion used  by  Robinson  [110]  and  by  Borgman  [21]  has  been 
modified  to  correspond  to  the  subscripts  values  acceptable 
to  FORTRAN  compilers.   The  Subroutine  NL0GN  approximates  the 
Fourier  integral  transform  pair  by  the  following  discrete 
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Fourier  series  approximation: 


LX 
n(n)  =  Af  Z   X(in)-exp  i{27T  (m-1)  (n-l)/LX}  (2.1a) 

m=l 

LX 
X(m)  =  At   E   n(n)-exp  i{-2TT  (m-l)  (n-l)/LX}  (2.1b) 

n=l 

where  the  time  sequence,  n (n) ,  is  given  at  LX  discrete  values 
of  time  which  are  partitioned  into  equal  intervals  of  At 
seconds  and  the  Fourier  spectral  sequence,  X(m) ,  is  given  at 
LX  discrete  values  of  frequency  which  are  partitioned  into 
equal  intervals  of  Af  Hertz.   In  order  for  the  inverse  trans- 
form to  be  exact,  the  following  relationship  between  the  time 
and  frequency  intervals  must  hold: 

At  •  Af  =  1/LX  (2.2) 

and  LX  must  be  equal  to  2  raised  to  an  integer  power. 

The  complex  Fourier  spectral  sequence  represented  by  Eq. 
(2.1b)  represents  a  two-sided  spectrum  with  one-half  of  the 
spectral  sequence  contained  in  the  complex  array  elements 

2  <  m  <  (LX+2)/2  (2.3) 

and  the  complex  conjugate  spectral  sequence  contained  in  the 
complex  array  elements 

(LX+2)/2  <  m  <  LX  (2.4) 

The  mean  value  of  the  time  sequence  is  represented  by  the 
first  array  element  (m=l).   The  time  sequence,  n  (n) ,  is  real 
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and,  therefore,  the  following  complex  conjugate  relationship 
for  the  complex  spectral  sequence  holds: 

X(LX+2-m)  =  X*(m)  (2.5) 

This  equality  may  be  demonstrated  by  substituting  m  = 

(LX+2-m)  into  Eq.  (2.1b)  to  obtain 

LX 
X(LX+2-m)  =  At   E   n(n)-exp  i{ -2tt  (LX-m+1)  (n-1) /LX}      (2.6) 
n=l 

The  exponential  argument  reduces  to  the  following: 

ARG  =  i{2Tr(m-l)  (n-l)/LX}-i{2TT(n-l)}  (2.7) 

which  gives  Eq.  (2.5)  as  a  result  of  the  following  identity: 

exp  i{-2Tr(n-l)}  =   1    ;    n=l ,  2 ,3. .  .  ,LX  (2.8) 

The  norming  constant,  LX,  is  associated  with  the  time 
sequence  in  the  Subroutine  NL0GN  (cf.  Robinson  [110],  pp. 
62-64)  vice  with  the  frequency  sequence  established  in  the 
standardized  notation  presented  in  Chapter  2.   Consequently, 
the  complex  Fourier  coefficients  for  the  linear  random  sea 
given  by  Eq.  (1.10)  were  initialized  by  the  following  complex 
expression: 

X(m)  =  F(m)-  (LX) 


=    (LX)    •    /SCmAa)    •    tt/T  '   •    exp    i{-il;(m)}  (2.9) 

BB  ^ 

where    ^Km)    is   the    random  phase    angle   generated  by   the   process 

described   in   Section   1   and  S(mAa)    is    the   value   of   the 

BB 
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Bretschneider  spectral  component  at  frequency,  m-ZTT'/Tp. 

R 

The    IBM  Subroutine    RANDU  was    used  to   generate   the    random 
number  sequence    required   for  the    computation  of  the    random 
phase    angles   uniformly   distributed   in   the    interval    {0,    Ztt}. 
A  description   of  the    simulation   of   a  linear  Gaussian   sea  using 
uniformly   distributed  phase    angles    is    given  by   Borgman    [20,23] 
and  by   Brown    [31] . 

3.      The   Bretschneider   Spectrum   and 
the   Phillips   Equilibrium   Spectrum 

Bretschneider    [28]    has    developed  the    following   spectral 
representation    for   a   linear   sea: 

2 
S(a)    =   ^-  exp{-0.675    [l]^}  (3.1) 

BB  a^  T  . 

where    a   and  T  are   parametric   constants.      This    one-sided  spec- 
trum  is    completely   defined  by   these   two  parameters.      In   order 
to    compare   Eq.     (3.1)    with   spectra   from   random  seas    obtained 
from  measured  hurricane    records,    it    is   necessary   to   transform 
these   two   constant   parameters    into   parameters   which    are 
measurable    from  the   hurricane    records.      This   may  be    accom- 
plished by    (1)    requiring   that   the   Bretschneider   spectrum   and 
the   measured  spectrum  have   equal   variance    and    (2)    that   the 
frequency   of  the  well-defined  peak   of  the   Bretschneider   spec- 
trum  agree    in   a  best   least-squares    sense  with   the    frequency 
of   the   often   ill-defined   "peak"   of  the   measured   spectrum.      As 
measured  spectra   frequently   differ   in   the   number   and  the   mag- 
nitudes   of   "peaks"  present    depending,    inter   alia,    on  the 
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degree   of  smoothing  used,    a  resort   to   a  best   least-squares 
agreement  between   the    frequencies    of   the  measured  and 
Bretschneider   "peaks"   is    required. 

The   variance   of   a  measured  spectrum  may  be    computed 
either    from  the   mean   square   of  the    recorded   realization, 

1        Tr/2      ^ 
E.   =   ^       /        n^(t)    dt  (3.2) 

^         ^R   -Tr/2 

where  Tn  is  the  length  of  the  measured  record  or,  equiva- 
lently,  from  the  integral  of  the  one-sided  spectral  density 
function 

oo 

E^  =  /  S(a)  da  (3.3) 

6  nn 

The  frequency  at  which  the  peak  of  the  Bretschneider  spectrum 

occurs  may  be  determined  by  equating  to  zero  the  derivative 

with  respect  to  the  angular  frequency  of  the  spectral  density 

function  and  evaluating  the  resulting  algebraic  equation  at 

the  frequency  at  which  the  peak  of  the  spectrum  occurs;  i.e., 

9  S(a) 

— 1^  =  0  ;   0  =  a  (3.4) 

3a  ° 

Equations    (3.3)    and    (3.4)    yield  two  equations   which   allow   the 
two   constant   parameters    given  by  Bretschneider  to  be   trans- 
formed into   two  measurable    constant  parameters,    i.e., 

{a,T}   -   {E^,a^}  (3.5) 

Integrating  the  Bretschneider  spectrum  is  most  easily  done  by 
first  making  the  following  change  of  variables: 
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T   a 


a 


(3.6a) 


V  =  exp{-0.675(^)^}  (3.6b) 

dv  =  ^  (|)^  exp{-0.675(^)^}  (3.6c) 

a 

Making  these  change  of  variables  and  integrating,  the  follow- 
ing expression  for  the  variance  of  the  spectrum  is  obtained: 

E.  =  -^4-ir   /  dv  =  -^4-r  exp{-0.675(i)4}  ? 
^    2.7(a)^  0      2.7(a)^  ^  o 

Ef  =   °^^  4  (3.7) 

^   2.7(0)"^ 

Differentiating  Eq.  (3.1)  with  respect  to  the  angular  fre- 
quency, a,    equating  the  result  to  zero,  and  evaluating  the 
resulting  expression  at  the  frequency  of  the  peak  in  the 
spectrum  yields  the  following: 

3  S(a)  _         2 

BB     r-5  .  2.7  rO    ^5,    as    ..._r  ^  ^nrrO    ^4- 


^^  ^  "^  ^S-J  >•  ^^  exp{-0.675(^)-};  a  =  a^ 


8a      "o    a       "o     (a„)"  ^o 

=0  (3.8) 

Equations  (3.7)  and  (3.8)  finally  determine  the  following  two 
constant  parameters: 

2.7  E, 


^  (a)^  (3.9a) 


2 
g 

(a)"^  =  2^  (a^)"^  (3.9b) 

Substituting  these  two  parameters  into  Eq.  (3.1)  results  in 
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the    following   form   for   the   Bretschneider  spectrum: 

S(a)    =   -— i   C^)^   exp{-1.25(^)^  (3.10) 

BB  ^o        "^  ^ 

The  best  least-squares  estimate  for  the  "peak  frequency,' 
a    ,   may  be  computed  by  the  application  of  the  linear  Taylor 
differential  correction  technique  (cf.  Marquardt  [86]  or 
McCalla  [87]).   In  this  technique,  a  mean  square  error 
between  the  measured  and  predicted  spectra  at  the  M  discrete 
frequencies  sampled  is  formed  as  follows: 

2    1   ^C  2 

e^  =  i-  Z   {S(m)-S(m)}'^  .,  ,,. 

c  m=l   nn   BB  l-^-J-J-J 

where  M  corresponds  to  the  cutoff  frequency  above  which  the 
measured  spectral  estimates  are  negligible. 

The  Bretschneider  spectral  distribution  function  is 
expanded  in  a  Taylor  series  about  the  best  estimate  of  the 
peak  frequency  in  terms  of  the  linear  differential  correction 
to  the  peak  frequency,  i.e., 

9    1   ^c  8  Sfm)  „  - 

e^  =  ^  E   {S(m)-S(m)-    .g^   (6a^) -0  (6a J  ^ ^         (3.12) 
c  m=l    nn    BB      ^o     °      ° 

Retaining  only   the    linear   corrections,    the   mean   square   error 
is    differentiated  with   respect   to   the    differential    correction 
to   the   peak   frequency,    <Sa    ,    and  the    result   equated  to   zero, 
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2  Mc  8    S(m)  d   S(m) 

^^=-2^     Z      (S(in)-S(m)-    -^^.(6a^)}.    -^^ 
^  *•      o^  c  in=l        nn        BB  o  o 


(3.13) 


which  yields    the    following   algebra  equation   for   6a    : 

M^-  8   S(m) 

Z      [{S(m)-S(m)}         .g^    ] 

6a      =   S^i 12 BB o_  ^3^^^^ 

°  M  8    S(m) 

in=l  o 

This  equation  may  be  iterated  for  the  (j+1)  correction  to 

the  i    estimate  for  &o      as  follows: 
-^  o 

M^  d   S(m) 

E   [{S(m)-s(^J'^m)}  jj^] 
m=l    nn   BB      '^'-^o^ 


(j+1) 


(<Sa  )^-'  ^^  =  ^^.j (3.15) 

°  M    9  S^J^(m) 


BB     ,2 


The  iteration  is  terminated  when  successive  corrections  are 
acceptably  small.   Table  3.1  demonstrates  the  results  of 
fitting  four  measured  spectra  from  Hurricane  Carla  with  a 
best  least -squares  estimate  for  the  peak  frequency.   The  fit 
and  number  of  iterations  to  convergence  are  the  poorest  for 
the  very  narrow  spectrum  represented  by  Record  No.  06886/2. 
The  number  of  spectral  estimates,  M  ,  used  from  the  measured 
spectra  was  equal  to  305  and  corresponded  to  a  cut-off  fre- 
quency equal  to  2.34  rad/sec  or  0.37  Hertz. for  Record  Nos . 
06885/1,  06886/1  and  /2 ,  and  to  2.92  rad/sec  or  0.47  Hz  for 
Record  No.  06887/1. 
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A  note  of  caution  is  in  order  with  regard  to  the  appli- 
cation to  design  of  the  Bretschneider  spectrum  expressed  as 
a  function  of  the  two  parameters  of  variance,  E-,  and  peak 
frequency,  a^.      Bretschneider  [27,28]  correlated  his  spectrum 
with  the  two  parameters  of  wind  speed  and  fetch  length  and 
there  is  no  satisfactory  guarantee  of  preventing  high  fre- 
quency instabilities  indicated  by  breaking  waves  in  the 
simulation  of  random  nonlinear  sea  realizations  when  the 
total  energy  content  and  peak  frequency  of  the  spectrum  are 
used.   In  the  absence  of  a  detailed  stability  analysis,  the 
equilibrium  spectrum  developed  by  Phillips  [99]  from  a  dimen- 
sional analysis  was  employed  to  limit  the  energy  content  in 
the  high  frequency  region.   Although  the  validity  of  the 
equilibrium  spectrum  has  been  questioned  [42,52,62,80,109], 
its  validity  for  application  to  engineering  design  appears 
to  be  appropriate. 

Phillips     (cf.  [100],  Ch.  4.5,  pp.  109-119)  reasoned 
by  similarity  conditions  that  the  spectral  components  in 
frequency  space  for  frequencies  greater  than  the  peak  fre- 
quency should  decay  according  to 

S(a)  =  3  •  g^/a^  ;  o    >   a  (3.16) 

nn 

Phillips    and  others    [99,100,104]    have    found  the    constant   6    to 
have    a  value    of  0.0117   ±    10%    for  a   expressed   in   radians   per 
second.       If  the    frequency   is   expressed   in  Hertz,    the    require- 
ment  that   the   energy   content   of  equivalent   differential    spec- 
tral   densities    remain   constant    implies    that 
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S(aJ-da  =  S(in/Tj^)-d(l/Tj^)  (3.17) 


(3.18) 


Tin                  nn 

which  yields 

6-g^.27T.(l/Tj^) 

B'-g^-(l/Tj^) 

[2TT-ni-(l/Tj^)]^ 

[m-d/Tj^)]^ 

so    that 

p.    _        B        _    . 

SI    X    10+10% 

p                    4         '  • 

(3.19) 


This  constant  for  the  equilibrium  spectrum  given  in  Hertz 
differs  from  that  given  by  Phillips  [100;  p.  114],  which 
appears  to  be  inverted.   The  nonlinearization  procedure 
includes  a  test  for  high  frequency  saturation  using  the  equi- 
librium spectrum  of  Phillips.   None  of  the  spectra  from  the 
nonlinear  simulations  of  the  measured  hurricane  spectra 
exceeded  the  equilibrium  spectrum. 

4.   Hurricane  Carla  Data  (September  8-10,  1961) 

Wave  Force  Project  II  (WPII)  recorded  the  instantaneous 
sea  surface  elevations  and  wave  pressure  forces  on  a  vertical 
piling  which  supported  an  instrumented  drilling  platform  at  a 
location  in  the  Gulf  of  Mexico  with  a  water  depth  of  approxi- 
mately 100  feet.   The  data  recorded  during  WPII  have  been 
made  available  to  the  public  through  the  National  Oceanographic 
Data  Center.   Thrasher  and  Aagaard  [122]  provide  details  of 
the  storms  recorded  and  other  information  relative  to  these 
data.   Hurricane  Carla  was  recorded  during  the  period 
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September  8-10,  1961.   The  continuous  records  which  are  avail- 
able from  this  period  of  WPII  contain  some  of  the  highest  and 
most  forceful  waves  recorded  during  WPII.   One  wave,  in 
particular,  had  an  average  trough  to  crest  wave  height  of 
almost  forty  feet.   These  data  contain  four  records  from 
Hurricane  Carla  with  continuous  sea  surface  elevations  and 
pressure  forces  recordings  each  covering  approximate  lengths 
of  time  from  eleven  to  fourteen  minutes.   The  data  are  digi- 
tized at  uneven  intervals  of  time  ranging  approximately  from 
0.12  to  0.24  seconds.   Calibration  information  and  equations 
for  transforming  the  data  into  engineering  units  are  given 
by  Blank  [16].   Table  3.2  contains  the  characteristics  of  the 
four  hurricane  records  used  in  the  comparison  with  simulated 
realizations.   The  data  were  redigitized  for  the  comparison 
analysis  by  linear  interpolation  after  applying  the  calibra- 
tion equations  given  by  Blank  [16]  in  order  to  obtain  records 
which  could  be  analyzed  by  the  Fast  Fourier  Transform  (FFT) 
algorithm.   No  time  shift  was  used  to  effect  coincidence 
between  the  sea  surface  realizations  recorded  by  the  wave 
staff  and  the  pressure  forces  recorded  on  the  instrumented 
pile  because  the  separation  distance  between  the  wave  staff 
and  center-line  of  the  pile  was  only  approximately  55  inches. 
Using  linear  wave  theory,  a  wave  with  a  period  of  ten  seconds 
traveling  in  an  ocean  of  uniform  depth  of  one  hundred  feet 
would  travel  the  straight-line  distance  between  the  center  of 
the  wave  staff  and  the  center  of  the  pile  in  approximately 
0.10  seconds.   Waves  incident  from  directions  other  than  alonj 
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Table    3.2 
Characteristics    of  Hurricane    Carla   Records 

Record   Length 
Record  No. Date Time (min.) LX  At 

06885/1  Sept.  9,    1961  2100  13.67  4096  0.20 

06886/1  Sept.  10,    1961  0000  15.01  4096  0.20 

06886/2  Sept.  10,    1961  0030  36.88  4096  0.20 

06887/1  Sept.  10,    1961  0300  11.00  4096  0.16 
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this    straight-line   would   require    less    travel   time    and,    there- 
fore,   the    time    shift  was   neglected.      The   effect   on  wave- 
induced  pressure    forces    of  neglecting   the    spatial    separation 
between   a  wave    staff  and   the    center-line   of  an   instrumented 
vertical   piling   is    demonstrated   in  Appendix  C  by   the    applica- 
tion  of  kinematics    from   linear  wave   theory   in   the  Morison 
equation. 

The  pressure    forces    for  each   individual   dynamometer  ele- 
vation  are    given   in   orthogonal    components  which   are    approxi- 
mately  aligned  with    the   platform  axes    (cf.    Blank    [16]    for 
dimensional    and   azimuthal    details).      The   pressure    forces 
recorded  by   the    dynamometer   at   the    55.33   foot   elevation  were 
used   in   the    comparison   analysis    since    this   elevation  was    con- 
tinuously  submerged   during   the   passage    of  waves    and  was    located 
near  enough    to   the    free    surface    to  have    recorded  the   effects 
of  the    smaller  high    frequency  waves   whose    influence    decays 
with   depth.      The    data   recorded   at    the    78.17    foot   elevation 
would  have   been  preferable    for   the   hydrodynamics    acting    at 
this   elevation  but    the    low   resolution   of  the   Y   force    component 

resulted   in  measured  values    changing   only  by   increments    of 

2 
approximately    30    lbs/ft      and  was    felt   to  be    less    accurate 

than   the    data   at    the    55.33    foot   elevation.      The    comparison 

between  measured   and  simulated  pressure    forces    involved  the 

resultant    forces    as    defined  by   Dean    and  Aagaard    [41].      The 

waves   were    assumed  to  be    collinear   and  to  propagate    in   a 

direction    9   =125°    relative    to   the   platform   coordinate    system. 
o 

The    resultant    force    at   time   nAt  was    obtained  by   the    following 


expression  : 


P(n)    =   sgn{Ae(n)}    •YFY(n)  +  FY(n) 
55  ^  ^ 
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(4.1) 


where  Fy(n)  and  FyCn)  are  the  orthogonal  pressure  force  com- 
ponents relative  to  the  platform  axes  digitized  at  time  nAt 
and   the   signum   function    is    determined  by   the    following: 


sgn{Ae(n)}    = 

where 

e(n)    =   ARCTAN 


+    for    |e(n)-e^|    <   90' 
-    for    |9(n)-e^|    >   90' 


FyCn) 


H^ 


(4.2) 


(4.3) 


The   effect   of   the    approximately    -4°    angular  misalignment  with 
respect   to   the   platform   coordinate    axes    of   the   dynamometer 
at   the    55.33    foot  elevation  was   not    included   in   the    determi- 
nation  of   the    direction   of   the    resultant    force. 


5.      Comparison   of   Sea   Surface   Realizations 

The   measured   spectra   from   four  Hurricane   Carla   records 
were   used   to   simulate   realizations    of   random  hurricane- 
generated  waves    correct    to   second  perturbation   order   for   two 
different   input    spectra.      The   smoothed  measured  spectra  of 
these    four   records    are   shown    in   Figs.    3.1,    3.2,    3.3,    and   3.4. 
Also   shown   on   these    figures    are    the   Bretschneider   spectra 
with   equal   variance    and  best   least-squares    fit    to   the   peak 
frequency   as    determined  by   the   procedure    discussed   in   Section    3. 
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The   measured  spectra  were   smoothed  by  block   averaging  over 
an    interval    of  nine    spectral   values.      This    type   of  smoothing 
filter  was    determined  by  Borgman    [21]    to  be    approximately 
equal   to   the    Gaussian   smoothing  which  hfe   used  to   investigate 
the    chi-squared   confidence    intervals    of  these    same   records 
(cf.    also   Bendat    and  Piersol    [6]).      Nine   spectral    components 
were   block    averaged  vice    the   eight    components    used  by  Borgman 
in   order   that   the    averaging  would  be   symmetrical.      The   block 
averaging  was   obtained  from  the   following: 

4 
S(m)    =      E        S' (m+n)/9  (5.1) 

nn         n=-4     nn 

where  S(m)  is  the  smoothed  measured  spectral  estimate  at  fre- 

nn 
quency  mAf  and  S'   is  the  raw  spectral  estimate.   The  discrete 

smoothed  spectral  estimates  have  been  connected  by  a  continu- 
ous curve  to  aid  in  reading  the  figures. 

Although  the  curve  which  connects  the  smoothed  measured 
spectral  estimates  appears  to  be  reasonably  smooth,  the 
erratic  nature  of  the  raw  spectral  estimates  caused  by  large 
differences  in  the  magnitudes  of  adjacent  estimates  and  by 
outlier  estimates  is  still  evident.   The  raw  spectral  esti- 
mates oscillate  rapidly  and  these  oscillations  are  especially 
prevalent  in  the  raw  spectra  near  the  best  least-squares  peak 
estimates  where  several  large  outlier  estimates  were  computed 
from  the  FFT  coefficients. 

The  variance  of  the  four  measured  records  varied  from 

2 
22.3-28.76  ft   and  the  best  least-squares  estimate  of  the  peak 
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frequency  varied   from  0.485-0.520    rad/sec.      The   root-mean- 
square   errors  were    computed   from  the    discrete   smoothed 
measured  spectral   estimates    and  the    discrete   Bretschneider 
spectral   values   by   the    following: 


>.=vti 


RMS  error  =  1/  ^  E   {S  (m) -S  (m)  l"-  (5.2) 

c  m=l   nn   BB 

where  M  is  the  total  number  of  spectral  estimates  used. 
Borgman  [21]  determined  that  the  magnitudes  of  the  raw  spec- 
tral estimates  became  negligible  after  three  hundred  and  five 
(305)  values  and  a  value  of  M  =305  was  selected  for  this  com- 
parison analysis.   The  RMS  errors  varied  between  5.29-12.23 
ft  /(rad/sec)  with  the  largest  error  occurring  for  the  rela- 
tively narrow  spectral  representation  for  Record  No.  06886/2. 
All  of  the  measured  spectra  demonstrated  several  "peaks"  and, 
except  for  Record  No.  06886/2,  there  are  two  "peaks"  of  nearly 
equal  amplitudes  very  near  the  best  least-squares  estimate  of 
maximum  peak. 

From  each  of  the  measured  and  Bretschneider  spectra 
representing  one  measured  record,  six  realizations  were  simu- 
lated.  Three  linear  and  three  nonlinear  realizations  for  both 
the  smoothed  measured  and  the  Bretschneider  spectra  were 
simulated  using  the  IBM  Scientific  Subroutine  Program  RANDU 
to  generate  the  random  numbers  required  to  compute  the  random 
phase  angles.   The  same  three  "seed"  numbers  required  to 
initialize  the  RANDU  random  number  generator  were  used  for 
each  of  the  four  records.   These  three  seeds  were  chosen 
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since  they  generally  demonstrated  the  following  three  differ- 
ent linear  measures  of  skewness:   (1)  negative,  (2)  positive, 
and  (3)  nearly  zero.   The  realizations  from  the  linear  simu- 
lation which  were  generated  using  this  method  would  have 
demonstrated  zero  skewness  measure  had  the  random  number  gener- 
ator realized  a  purely  Gaussian  process. 

The  cumulative  probability  distributions  from  each  of 
the  four  records  showing  all  six  simulated  realizations  in 
addition  to  the  measured  realization  are  shown  in  Figs.  3.5, 
3.6,  3.7,  and  3.8.   Realizations  from  the  unsmoothed  measured 
spectrum  would  have  demonstrated  exactly  the  same  distribution 
as  the  measured  data  represented  by  the  open  circles  if  the 
identical  measured  phase  angles  had  been  used  as  a  result  of 
the  exact  inverse  relationship  in  the  FFT  algorithm.   The 
variations  between  the  distributions  from  the  measured  spectrum 
demonstrate   the  importance  of  these  phase  angles  in  deter- 
mining the  distribution  of  any  realization.   The  straight 
line  represents  a  normal  distribution  having  the  identical 
mean  and  standard  deviation  as  the  realizations. 

All  realizations  demonstrate  close  agreement  to  a 
Gaussian  process  between  two  normalized  standard  deviations. 
The  ordinate  values  have  been  normalized  by  the  standard  devi- 
ation of  the  simulated  or  measured  record;  i.e.. 


t,(i)  =  jniil i  =  l,2,...,31    (5.3) 
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where  nCi)  is  the  value  o£  the  i^   discrete  cumulative  inter- 
val, o^   is  the  standard  deviation  of  the  measured  or  simu- 
lated realization;  and 


where  n    is  the  maximum  value  of  the  sea  surface  realiza- 
max 

tion  and  n  .   is  the  minimum  value.   The  mean  of  the  measured 
min 

or  simulated  record  was  subtracted  from  the  value  of  each 
realization  and  the  variance  was  computed  by 


rrXTij  Jl  {n(n)-n}2  (5.5) 


where  n  is  the  mean  of  the  realization  determined  by 

T   LX 
^  =  1   Z  n(n)  (5.6) 

^^  n=l 

and  LX  is  the  number  of  time  values  in  the  random  sequence. 

The  largest  values  of  the  sea  surface,  n ,  from  the  nonlinear 

simulations  are  seen  from  these  figures  to  be  approximately 

0.5-1.5  standard  deviations  greater  than  the  corresponding 

largest  realizations  from  the  linear  simulations. 

Figures  3.9,  3.10,  3.11,  and  3.12  demonstrate  the  second 

order  spectra,  -^S(a),  computed  for  both  the  smoothed  measured 

^  nn 
and  Bretschneider  spectra.   Theoretically,  these  spectra  may 

be  added  linearly  to  the  first  order  spectra  as  a  consequence 
of  the  Gaussian  assumption  and  the  resulting  total  energy,  E^, 
would  be  statistically  identical  for  all  nonlinear  realiza- 
tions from  the  same  initial  linear  spectra.   The  final  value 
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of  the    total    energy    for  each   realization   from  the    same 
initial   spectrum  varies    as    a   result    of   the   nonvanishing 
three-product   second  order  statistical   moment  between   the 
first   and  second  perturbation   orders    for  the    random  sea  sur- 
face,   i.e.,    the    linear   first  perturbation   order  solution   is 
not   strictly   Gaussian   and 

£{in(t)-2n(t+T)+2n(t)-^n(t+T)}  i  0  (5.7) 

The  discrete  second  order  spectra  may  be  computed  from 
the  continuous  first  order  spectra  given  by  Eq.  (6.8)  in 
Chapter  2,  Section  6,  by  replacing  the  integral  with  a  dis- 
crete sum:  i.e.. 


2S(m)  =  2   E   S(£-m)-S(m) •H^[(£-m),m]-Aa  (5.8) 

nn      ^=1  nn     nn 

where 

ha   =   p-  (5.9) 

Substituting  the  FFT  expressions  for  the  spectra  yields 

^c 
,S(m)  =  — i^  I      X(Ji-m)-X(m)-X*(il-m)-X*(m)  (5.10) 

nn     (LX)-^TT  £=1 

where  M  is  the  maximum  number  of  discrete  FFT  coefficients 
required  to  represent  the  energy  spectrum. 

The  second  order  spectra  from  both  the  smoothed  measured 
and  Bretschneider  spectra  consistently  indicate  a  peak  at 
approximately  2a  .   The  peak  spectral  values  from  the  smoothed 
measured  spectra  are  consistently  larger  in  magnitude  than 
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the  spectral  values  of  the  peak  estimates  from  the  Bret  Schneider 
spectra.   While  the  Bret Schneider  spectra  contain  only  one 
second  order  peak,  the  measured  spectra  generally  contain  more 
than  one.   These  second  order  peaks  reflect  the  nonlinearities 
of  the  self -interactions .   The  contributions  to  the  nonlinear 
realization  from  these  self-interactions  indicate  the  impor- 
tance of  the  outlier  spectral  estimates  which  are  capable  of 
generating  large  self-interaction  values.   Both  the  smoothed 
measured  spectra  and  the  Bretschneider  spectra  contain  smoothly 
varying  spectral  estimates  and  the  possibility  of  large  non- 
linear self -interaction  contributions  has  been  eliminated  as 
a  result  of  this  smoothness. 

Ensembles  from  three  linear  and  nonlinear  realizations 
from  the  smoothed  measured  spectra  are  shown  in  Figs.  3.13, 
3.14,  3.15,  and  3.16.   Ensembles  from  three  linear  and  non- 
linear realizations  from  the  Bretschneider  spectra  are  shown 
in  Figs.  3.17,  3.18,  3.19,  and  3.20.   The  time  intervals  of 
the  measured  realizations  shown  at  the  top   of  these  ensembles 
were  selected  because  the  largest  waves  measured  were  recorded 
during  these  intervals.   The  simulated  ensembles  generally  did 
not  contain  the  largest  waves  in  the  simulated  realization  and 
the  figures  of  cumulative  distributions  must  be  examined  in 
order  to  determine  the  extreme  values  which  were  realized  for 
any  particular  simulation. 
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Figure  3.13.   Ensemble  Comparison  Between  Measured  Realiza- 
tion and  Linear  and  Nonlinear  Realizations 
Simulated  from  Smoothed  Measured  Spectrum  from 
Record  No.  06885/1. 
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Figure    3.14.      Ensemble   Comparison   Between  Measured   Realiza- 
tion  and   Linear   and  Nonlinear   Realizations 
Simulated   from  Smoothed  Measured   Spectrum   from 
Record  No.    06886/1. 
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Figure  3.15.   Ensemble  Comparison  Between  Measured  Realiza- 
tion and  Linear  and  Nonlinear  Realizations 
Simulated  from  Smoothed  Measured  Spectrum  from 
Record  No.  06886/2. 
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Figure    3.16.      Ensemble    Comparison  Between  Measured  Realiza- 
tion  and  Linear   and  Nonlinear  Realizations 
Simulated   from  Smoothed  Measured   Spectrum   from 
Record  No.    06887/1. 
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Figure  3.17.   Ensemble  Comparison  Between  Measured  Realiza- 
tion and  Linear  and  Nonlinear  Realizations 
Simulated  from  Bretschneider  Spectrum  from 
Record  No.  06885/1. 


101 


12 

|— 

Seed    No    7071 

— , 

- 

A 

h  ~ 

6 

"  A     /^ 

A         1 

1  \  - 

0 

■    \   / 

A              1       ^^           /^                /         \ 

(\       i 

\- 

/    \  /' 

/        ^      A 

\ 

-/    \  /' 

W            \    / 

\- 

-6 
.1? 

i     1 

1              1              1              1              1              1              1 

1      1      1 

^ 

533B 


543.8 


553.8 


563.8  573.8 

Time    (sec) 


583.8 


593.8 


6038 


Figure    3.18.      Ensemble   Comparison   Between  Measured   Realiza- 
tion  and  Linear   and  Nonlinear   Realizations 
Simulated   from  Bretschneider   Spectrum   from 
Record  No.    06886/1. 
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Figure  3.19.   Ensemble  Comparison  Between  Measured  Realiza- 
tion and  Linear  and  Nonlinear  Realizations 
Simulated  from  Bretschneider  Spectrum  from 
Record  No.  06886/2. 
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Figure  3.20.   Ensemble  Comparison  Between  Measured  Realiza- 
tion and  Linear  and  Nonlinear  Realizations 
Simulated  from  Bretschneider  Spectrum  from 
Record  No.  06887/1. 
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6.   Digital  Linear  Filter 

The  digital  linear  filter  technique  is  a  method  by  which 
certain  wave  fields  may  be  computed  by  the  convolution  of  an 
impulse  response  function  with  a  time  sequence  of  the  sea  sur- 
face.  Reid  [106]  developed  the  technique  to  compute  kinematic 
fields  for  predicting  wave  forces  on  piles.   Grooves  [53]  has 
tabulated  a  number  of  frequency  response  functions  using 
linear  wave  theory  which  will  compute  several  wave  fields  by 
this  general  convolution  method.   Wheeler  [129]  applied  the 
technique  developed  by  Reid  [106]  to  hurricane-generated  wave 
data  recorded  by  Wave  Project  II  during  Hurricane  Carla  and 
computed  drag  and  inertia  force  coefficients  which  varied  with 
the  vertical  depth  of  pressure  force  correlation.   A  succinct 
and  lucid  description  of  the  technique  which  is  applicable  to 
wave  force  predictions  has  been  given  by  Borgman  [22]. 

As  applied  in  the  prediction  of  wave  forces  on  piles,  the 
digital  linear  filter  technique  assumes  that  the  kinematic 
wave  fields  may  be  computed  from  linear  wave  theory  and  that 
the  sea  surface  realizations  represent  an  irregular  long- 
crested  wave  field  which  is  the  result  of  the  superposition 
of  linear  waves.   The  application  which  is  derived  in  this 
section  departs  from  the  descriptions  given  by  Reid  [106]  and 
by  Borgman  [22]  in  the  following  areas: 

(1)  a  "stretched"  vertical  coordinate  is  to  be 
utilized. 
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(2)  the  horizontal  acceleration  field  is  to  be  com- 
puted by  numerical  differentiation  of  the  hori- 
zontal velocity  field  vice  by  the  convolution 
method. 
The  stretching  of  the  vertical  coordinate  in  the  argument 
of  the  vertical  coordinate  dependent  function  in  the  expres- 
sion for  the  linear  velocity  potential  has  been  employed  with 
some  success  [67,129]  to  ensure  that  this  fimction  does  not 
become  excessively  large  in  magnitude  whenever  kinematic 
fields  are  evaluated  for  wave  crest  phase  values  at  eleva- 
tions above  the  still  water  level.   For  kinematic  predictions 
required  at  a  vertical  elevation,  s,  above  the  bottom,  the 
stretching  function  results  in  the  computations  actually 
being  made  at  a  stretched  elevation,  s',  above  the  bottom; 
i.e.  , 

s'(a,s)  =  a(h,t,x)-s  (6.1) 

where    the    vertical    stretching   function,    a(h,t,x),    is    given  by 
the    following: 

°'^h'^'^^  =  h.nU,xJ  C6.2) 

A  linear  velocity  potential  when  modified  by  this  stretch- 
ing term  becomes 

Substituting  Eq.  (6.3)  into  the  equation  of  continuity  yields 
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where  the  terms  0{*}  are  of  higher  order  in  wave  slope,  Hk/2, 
than  the  two  derivatives  on  the  left  side.   Eq.  (6.4)  demon- 
strates that  the  linear  velocity  potential  modified  by  a 
stretched  vertical  coordinate  no  longer  satisfies  exactly  the 
equation  of  continuity;  however,  the  modified  velocity  poten- 
tial does  satisfy  the  equation  of  continuity  to  first  order 
in  wave  slope  and  is,  therefore,  consistent  with  the  pertur- 
bation expansion  which  is  assumed  to  be  valid  for  linear  wave 
theory.   Additional  expressions  which  may  be  utilized  to 
extend  the  wave  pressure  field  computed  by  linear  wave  theory 
above  the  still  water  level  are  given  by  Bergman  [18], 
Chakrabarti  [33],  and  Nath  and  Felix  [93]. 

The  following  brief  derivation  of  the  linear  filter  tech- 
nique which  is  applicable  to  computing  horizontal  kinematic 
flow  fields  in  a  stretched  vertical  coordinate  follows  closely 
the  derivations  given  by  Reid  [106]  and  by  Bergman  [22]  and 
is  included  only  for  convenience. 

From  the  following  aperiodic  linear  velocity  potential 

$Cx,s',t)  =  -i  /  F(a)-  f  gg^};{j;gj^  exp  i{at-kx}da      (6.5) 

and  linear  dispersion  equation 

a^  =  gk  tanh{kh}  (6.6) 

With  k(-a)=-k(a),    the    following  horizontal   kinematic   fields 
and  scalar   sea   surface    realization   for   a   fixed  horizontal 
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spatial   value    for  x    (x=0,    say)    may  be   easily   computed   correct 
to   first    order   in  wave    slope: 

u(s',t)    =      f   F(a)  a^?^};{j^gj^    exp    i{at}d0  (6.7a) 


u^(s',t)    =    i      /   F(a)a2    sinMkh}^   ^^P    i{at}da  (6.7b) 


ri(t)    =      /   F(a)    exp   i{at}da  (6.7c) 

-co 

The   spectral   density   o£  the   sea  surface,    F(a) ,    is    obtained  by 
the    following   Fourier   inversion: 

CO 

F(a)    =  ^     /  n(t)    exp   i{-at}dt  (6.8) 

The    spectral    response    functions    for   the  horizontal   kinematic 
fields    are    seen  by   inspection   of  Eqs .     (6.7a,b)    to  be    the 
following: 

^     r^         t^  ^     COSh{ks'}  /■£     r>     -v 

^u^'^^'^^    =   ^    sinhlkhl  ^^'^^^ 

r         fn     oM     -     rr^     COSh{ks'}  ..     q,   >. 

%^^'^  ■>    -   ^   sinhlkhl  f^-^b^ 

which  are  symmetric  and  antisymmetric  functions,  respectively, 
of  the  frequency,  a,  provided  that  k(-a)=-k(a). 

These  spectral  response  functions  may  be  expanded  in 
generalized  Fourier  series  (cf.  Lighthill  [75],  Papoulis  [96], 
or  Titchmarsh  [126])  where  the  interval  of  periodicity  must  be 
greater  than  or  equal  to  the  maximum  frequency  of  measurable 
energy  in  the  spectrum  of  the  sea  surface,  F(a).   Reid  [106] 
did  not  employ  a  stretched  vertical  coordinate  and  truncated 
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the   frequency  response   functions    at   a  cut-off  frequency  which 
was    less    than   the    frequency  of  periodicity.      IVheeler    [129] 
used   a   stretching   fionction   and   let   the   value   of  the    cut-off 
frequency  equal    the   periodic   interval.      The    cut-off   frequency 
must  be    less    than   or  equal    to   the   Nyquist    frequency  which    is 
determined  by   the    spacing   used   in   digitizing   the   measured 
record;    i.e., 

^cik  C6.10) 

where    At   is    the    digitizing   interval    and  a      is    the    cut-off 

angular   frequency. 

Expanding   the    frequency   response    functions    given  by   Eq. 

(6.9a,b)    in   generalized   finite    Fourier   series  which   are 

periodic   in    frequency   space    over   the    interval    of  periodicity, 

{2tt/t    }      ,    the    following   series    result: 

M 
G^ia,s')    =      Z        C^  exp   i{-amT^}  (6.11a) 

m=-M 

M 
^u    (^^'^'^    =    12        C^  exp   i{-amT    }  (6.11b) 

t  m=-M  ° 

where   M   is    the   maximum  number   of  coefficients    required  to   fit 
the    series    to   the    frequency   response    function,    {Ztt/t    }        is 
the    interval   of  periodicity,    and   the    reality   requirements    for 
symmetric   and   antisymmetric  expansions    are   given   by   the 
following   relationships: 

^-m  ^   ^m  (symmetric)  (6.12a) 

^-m  ^  "^m      (antisymmetric)  (6.12b) 
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The  generalized  Fourier  transforms  of  Eqs .  (6.11a,b) 
(cf.  Lighthill  [75],  Papoulis  [96],  or  Titchmarsh  [126])  are 
given  by  the  following  impulse  response  functions: 

00 

g    (t,s')  =  /  g  CcJ,s')  exp  i{aT}da  (6.13a) 

g      (t,s')  =  /  g^  CcJ,s')  exp  i{aT}da  (6.13b) 

t  -00      t 

Substitution  of  the  Fourier  series  approximations  for  the 

frequency  response  functions  given  by  Eqs.  (6.11a,b)  into  Eqs. 

(6.13a,b)  yields  the  following: 

M       «> 
g    (t,s')  =   Z    C   /  exp  i{a(T-mT^)}da  (6.14a) 

m=-M     -<=° 

M 
g      (t,s')  =i  E   C   /  exp  i{a(T-mT^)}da  (6.14b) 

^t         m=-M  ^   -^ 

The  integrals  in  Eqs.  (6.14a,b)  are  seen  to  be  equal  to  the 
Dirac  delta  distribution  in  the  time  domain  by  their  similarity 
to  Eq.  (2.12),  Chapter  2,  Section  2,  which  defines  the  shift- 
ing property  in  the  spatial  domain.   Therefore,  the  impulse 
response  functions  are  seen  to  be  a  comb  of  Dirac  delta  dis- 
tributions in  the  time  domain  and  are  proportional  to  the 
Fourier  series  coefficients  which  approximate  the  continuous 
frequency  response  functions,  i.e., 

M 
g^(T,s')  =  2tt   E   V(T-mT^)  (6.15a) 

m=-M 


g      (t,s')  =  27ri   E   C^6(T-mTQ)  (6.15b) 

t  m=-M 
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The   expressions    for   the   Fourier   transform  o£   the   sea 
surface   spectral    density,    F(a),    given  by  Eq.     (6.8)    and  the 
expressions    for   the    frequency   response    functions    given  by 
Eqs.    (6.9a,b)    may  be    substituted  into   the   equations    for   the 
aperiodic  horizontal   kinematic   flow    fields    given  by   Eqs. 
(6.7a,b)    to  yield  the    following: 

00  00 

u(s',t)    =    /   G^(a,s')-exp    i{at}-    ^7     ^  ri(T)exp   i{-aT}dTda 

-0°  -<»  (6.16a) 

00  00 

u.  (s',t)    =    i    /   G,,    (a,s')-exp    i{at}-    7—     /  n(T)exp   il-atldida 
^  -00     "t  ^"^    -co  (6.16b) 

Inverting   the   order  of   integration   and  rearranging   terms 

results    in   the    following: 

u(s',t)    =    /   ti(t)-{^     /   G^(a,s')exp   i{    a(t-T)}da}dT  (6.17a) 

-00  -00 

00  00 

u^(s',t)    =    i    /   n(T)-{27     /   G^    (a,s')exp    i{a(t-T)}da}dT    (6.17b) 
-  00  -00        t 

Comparing  the  expressions  in  brackets  {•}  in  Eqs.  (6.17a,b) 
with  the  expressions  for  the  impulse  response  functions  given 
by  Eqs.  (6.13a,b),  it  may  be  seen  that  these  bracketed  terms 
are  the  impulse  response  fimctions  shifted  by  (t-x).   Replac- 
ing T  with  (t-r)  in  the  equivalent  expression  for  the  impulse 
response  function  given  by  Eqs.  (6.15a,b)  and  substituting 
for  the  bracketed  terms  in  Eqs.  (6.17a,b),  the  following  are 
obtained: 

M 
u(s',t)  =   Z   C   /  n(T)-6[t-(T+mT^)]dT  (6.18a) 

m=-M  "^  -"  ° 

M       " 
u^(s',t)  =  -E    C'   /  n(T)-5[t-(T+mT^)]dT  (6.18b) 

t         ni=-M  "^  -co 
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The  following  change  of  variables  substitution 

T'  =  t-T  (6.19) 

results  in  the  final  convolution  expressions  for  the  horizontal 
kinematic  fields: 

M  00 

u(s',t)=   E   C   /  n(t-T')[6(T'-mT  )  +  6(T'+mTQ)]dT'      (6.20a) 
m=0     -oo 

M 
u^(s',t)  =   Z   C;   /  n(t-T')[5(T'+inTQ)-6(T'-inTQ)]dT'    (6.20b) 
in=l     -oo 

The  convolutions  given  by  Eqs.  (6.20a,b)  are  the  weighted 
sums  of  products  of  the  sea  surface  sequence  with  the  Fourier 
coefficients  of  the  spectral  response  functions. 

These  spectral  response  functions  are  a  function  of  time 
through  the  stretching  fionction.   Consequently,  the  spectral 
response  function  for  the  horizontal  velocity  was  partitioned 
into  ten  even  intervals  of  constant  elevation,  s',  given  by 


J  =  £-s'/10  a    =    0,1,... ,10       (6.21) 


and  the  instantaneous  horizontal  velocity  was  computed  at  the 
eleven  contours  of  constant  elevation,  s'   from  Eq.  (6.20a). 
The  instantaneous  horizontal  velocity  at  the  desired  elevation, 
s,  was  then  computed  by  linear  interpolation  between  contours 
of  constant  elevation,  s' 

The  previous  applications  of  the  vertical  "stretching" 
function  involved  predicting  pressure  forces  only  for  maximum 
peak  forces  [129]  or  only  for  forces  recorded  during  the  crest 
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portion  of  a  wave  [67].   The  advantages  of  extending  the 
kinematics  of  linear  theory  to  vertical  elevations  above  the 
still  water  level  by  Eq,  (6.2)  are  not  achieved,  however, 
without  a  corresponding  deleterious  eff*ect  of  amplifying 
kinematic  field  predictions  in  the  trough  region  of  a  wave 
compared  to  predictions  from  unstretched  linear  theory. 
Since  a  continuous  time  sequence  of  predicted  pressure  forces 
are  required  for  simulation,  this  overprediction  by  the 
stretching  function  in  the  trough  region  of  waves  can  be  as 
important  as  the  overprediction  in  the  crest  region  by  the 
unmodified  linear  theory.   A  demonstration  and  discussion  of 
this  effect  may  be  found  in  Section  7  of  this  chapter. 

In  the  application  to  the  WPII  data,  the  cut-off  fre- 
quency of  0.3125  Hz  and  the  filter  spacing  of  1.60  seconds 
used  by  Wheeler  [129]  were  adopted.   The  instantaneous  hori- 
zontal velocity  for  a  desired  elevation  was  computed  from  the 
convolution  given  by  Eq.  (6.20a)  using  the  weighting  coeffi- 
cients for  contours  of  constant  elevation  shown  in  Table  3.3 
and  by  linear  interpolation  between  these  contours  of  constant 
elevation.   The  instantaneous  horizontal  acceleration  was 
computed  by  numerically  differentiating  the  horizontal  veloc- 
ity at  contours  of  constant  elevation,  s^,  according  to  the 
following  formula: 

u[s^,n+l]-u[s^,n-l] 


•H  rt 
J"  C 

O-H 
Jh   O 

o  o 
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and  the  value  of  the  horizontal  velocity  at  the  desired  ele- 
vation, s  ,  was  determined  by  numerical  interpolation  between 
contours  of  constant  elevation.   The  resultant  horizontal 
pressure  on  a  vertical  pile  at  elevatioh,  s,  was  computed 
from  the  Morison  equation  utilizing  these  two  kinematic 
fields  and  the  modified  Dean  and  Aagaard  resultant  force 
coefficients  from  Table  3.4  (cf.  [67]).       The  Reynolds 
number  was  computed  using  a  value  for  the  kinematic  viscosity 
equal  to  1.4x10    ft  /sec  and  for  a  piling  diameter  equal  to 
3.71  ft.   The  density  of  the  sea  water  was  assumed  to  be  2.0 
slugs/ft^. 


Table  3.4 

Drag  and  Modified  Inertia  Force  Coefficients 

for  Resultant  Pressure  Forces 

(from  Dean  and  Aagaard  [41]) 


Reynolds 
(xlO'^) 

No. 

Re    <    3 

3   <    Re    <    10 

10    <   Re 

^D 

1.34 
1.18 

0.98 
1.18 

0.92 
1.18 

7.   Comparison  of  Pressure  Forces 

To  reiterate,  the  main  purpose  of  the  nonlinear  random 
sea  simulation  is  to  provide  a  random  forcing  function  to  be 
employed  in  the  stochastic  analysis  of  deep  water  pile- 
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supported  platforms.   The  principal  comparisons  between  the 
random  sea  simulations  and  the  measured  realizations  which 
were  previously  discussed  involved  only  the  statistics  of  the 
random  waves.   Borgman  [17]  has  developed  a  theory  for  deter- 
mining wave  forces  on  a  pile  from  a  narrow  band  spectrum  of 
random  waves.   Wells  [128]  offered  constructive  comments  rela- 
tive to  the  second  order  nonlinearities  omitted  in  the  theory. 
Pierson  and  Holmes  [102]  later  developed  a  probability  den- 
sity function  for  random  wave  forces  using  Gaussian  assump- 
tions and  compared  both  the  probability  densities  and  cumula- 
tive probability  distributions  of  their  model  with  measured 
force  data.   The  advantages  and  disadvantages  of  the  model 
were  discussed  by  several  writers  [40,90,107]  in  addition  to 
the  closure  by  the  authors  [103]. 

No  attempt  has  been  made  in  this  study  to  derive  the 
probability  density  function  of  the  pressure  force  distribu- 
tions obtained  by  filtering  the  nonlinear  random  realizations 
as  this  is  a  most  complex  problem  due  to  the  nonlinearities 
involved.   One  simulated  realization  from  each  of  the  four 
measured  hurricane-generated  continuous  records,  however, 
was  selected  as  a  random  input  for  the  digital  linear  filter 
technique  in  order  to  compute  pressure  forces  for  comparison 
and  demonstration  of  the  force  realizations  which  may  be 
expected  in  applications.   The  same  seed  number  for  the  random 
number  generator  was  used  for  each  simulated  realization  since 
the  sea  surface  realizations  computed  from  the  same  seed 
number  represented  varying  degrees  of  fits  to  distributions 
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o£  the  measured  sea  surface  realizations  from  each  measured 
record  and,  therefore,  did  not  always  give  the  best  agree- 
ment with  the  measured  realizations. 

The  stretching  function,  a(h,t,x),  introduced  in  Section  6 
is  inversely  proportional  to  the  instantaneous  value  of  the 
sea  surface,  n.   This  fact  introduces  a  negative  bias  into 
the  computation  of  horizontal  kinematics  by  the  digital  linear 
filter  technique  which  results  in  the  overprediction  of  kine- 
matics in  the  trough  region  of  the  waves  compared  to  pre- 
dicted kinematic  values  from  unstretched  linear  theory. 

Figures  3.21  and  3.22  illustrate  this  effect  for  two 
strictly  periodic  waves.   The  nonlinear  wave  shown  by  the 
dotted  line  is  a  dimensional  example  of  the  Dean  Stream  Func- 
tion CASE  7-C  and  represents  a  nondispersive  progressive  wave 
having  a  wave  height  which  is  three -fourths  of  the  breaking 
wave  height  for  a  water  depth  of  99  feet  and  a  wave  period 
of  10  seconds.   The  solid  line  represents  a  linear  wave  with 
the  equivalent  wave  period  and  total  energy  as  the  dimensional 
form  of  the  Stream  Function  CASE  7-C.   The  horizontal  veloci- 
ties were  computed  using  the  filter  coefficients  given  in 
Table  3.3  and  the  horizontal  accelerations  were  computed  by 
numerical  differentiation  using  Eq.  (6.22).   The  pressure 
forces  were  computed  by  the  Morison  equation  using  the  modi- 
fied Dean  and  Aagaard  coefficients  from  Table  3.4.   The  values 
of  the  kinematic  viscosity,  sea  water  density,  and  piling 
diameter  required  to  compute  the  Reynolds  number  and  pressure 
forces  are  given  in  Section  6.   The  negative  bias  in  the 
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pressure  force  is  noticeable  in  both  examples  with  the  more 
dramatic  effect  demonstrated  at  the  55.3  foot  elevation.   The 
effect  of  the  skewed  Stream  Function  profile  appears  to  be 
lower  velocity  values  compared  with  the  value  computed  by- 
filtering  the  unskewed  linear  profile  in  the  trough  regions 
of  the  two  profiles. 

Figures  3.23,  3.24,  3.25,  and  3.26  each  contain  four 
pressure  force  spectra  from  the  dynamometer  located  at  the 
55.3   foot  elevation.   The  four  spectra  represent  the  follow- 
ing:  (1)  measure  pressure  force  spectra;  (2)  pressure  force 
spectrum  computed  from  filtering  a  nonlinear  simulated  reali- 
zation synthesized  from  the  measured  amplitude  spectrum; 
(3)  pressure  force  spectrum  computed  from  filtering  a  non- 
linear simulated  realization  synthesized  from  the  Bretschneider 
spectrum;  and  (4)  pressure  force  spectrum  computed  from  filter- 
ing the  measured  sea  surface  realization. 

Table  3.5  demonstrates  the  comparisons  of  mean  square 
pressures  determined  by  these  four  methods  for  each  of  the 
four  measured  records.   None  of  the  predicted  mean  square 
pressures  exceeded  the  measured  values.   The  mean  square 
pressures  computed  from  the  measured  realizations  ,l/Pp  ,were 
consistently  greater  than  the  other  predicted  values.   The 
mean  square  pressures  predicted  from  realizations  simulated 
from  the  Bretschneider  spectrum,\/Ppg,  were  consistently 
greater  than  the  values  from  realizations  simulated  from  the 
smoothed  measured  spectrum,l/Ppj.;  however,  these  differences 
were  not  great. 
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In  general,  the  shape  of  the  spectra  computed  from  filter- 
ing the  measured  realization  approximated  very  closely  that  of 
the  measured  sea  realization  and     force  spectra.   This 
similarity  was  noticed  by  Wiegel  and  le*d  Borgman  [18]  to 
develop  his  linear  spectral  model.   It  is  possible  that  an 
adjustment  of  the  force  coefficients  could  bring  these  predicted 
pressure  force  spectra  into  near  perfect  agreement  with  the 
measured  pressure  force  spectra.   In  any  event,  it  is  evident 
that  the  force  spectra  computed  from  the  measured  realizations 
agree  very  well  with  the  measured  force  spectra. 

Figures  3.27,  3.28,  3.29,  and  3.30  demonstrate  the  pres- 
sure force  distributions  from  the  following  four  types  of 
measured  and  simulated  realizations:   (1)  measured;  (2)  pre- 
dicted from  filtering  the  measured  sea  surface  realization; 
(3)  predicted  from  filtering  both  the  linear  and  nonlinear 
simulated  realization  synthesized  from  the  smoothed  measured 
spectrum;  and  (4)  predicted  from  filtering  both  the  linear 
and  nonlinear  simulated  realizations  synthesized  from  the 
Bretschneider  spectrum.   Again,  the  same  seed  number  for  the 
random  number  generator  was  used  for  all  simulations.   As  was 
noted  with  the  distributions  of  the  random  sea  surface  reali- 
zations, the  measured  pressure  force  data  and  the  pressure 
forces  computed  from  simulated  realizations  agree  very  well 
between  two  normalized  standard  deviations.   The  ordinate 
values  in  these  figures  have  been  normalized  by  the  standard 
deviation  of  the  realization  by  the  following  expression: 
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P55  (i)  =    ^^  "^'\ —      ;   i  =  l,2,...,31        (7.1) 

where  Pr^  (i)  is  the  normalized  value  of  the  i    discrete 
cumulative  pressure  force  interval,  o^   is  the  variance  of  the 
measured  or  simulated  pressure  force  realization  at  the  55.3 
foot  dynamometer  elevation,  and 

^^SS   -    [fsS^max-Css'minl/^''  (^-2) 

■where    (Pcc^j^ax   ^^    ^^^  maximum  pressure    force    and    (Prr)    •      is 
the   minimum  pressure    force   measured  or  predicted   at   the    55.3 
foot   dynamometer  elevation.      The   mean   of  the   pressure    force 
time    sequence  was    subtracted   from   the   value   of  each    sequence 
and   the   variance   was    determined   from   the   following: 

LX 

where  Prr  is  the  mean  of  the  sequence  determined  from 

T   LX 
F   =  f^  Z   P(n)  (7.4) 

^^      n=l   55 

and  LX  is  the  total  number  of  values  in  the  random  sequence. 

Table  3.6  contains  the  statistical  measures  from  each  of 
the  four  pressure  force  time  sequences  obtained  from  the  four 
measured  records.   The  skewness  of  each  of  the  four  measured 
records  was  less  than  0.2.   Record  No.  06885/1  demonstrates 
a  negative  skewness  measure  for  both  the  sea  surface  and 
pressure  force  realizations.   Three  of  the  measured  pressure 
force  records  demonstrate  kurtosis  measures  greater  than  1.40. 
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The  low  negative  kurtosis  measure  came  from  the  relatively 
narrow  band  spectrum  representing  Record  No.  06886/2. 

Both  the  skewness  and  kurtosis  measures  of  the  pressure 
forces  computed  by  filtering  the  measured  sea  surface  reali- 
zations agree  favorably  with  these  same  measures  from  the 
measured  pressure  force  time  series.   This  agreement  is 
another  indication  of  the  accuracy  of  the  digital  linear 
filter  technique  in  predicting  pressure  forces  provided  the 
appropriate  nonlinear  random  realization  is  filtered. 

A  comparison  of  the  predictions  from  simulated  realiza- 
tions indicates  that  the  predictions  from  the  linear  simula- 
tions demonstrate  consistently  negative  skewness  measures 
whose  magnitudes  were  relatively  unaffected  by  the  shape  of 
the  spectrum  used  for  the  simulation.   None  of  the  nonlinear 
predictions  demonstrate  negative  skewness  measures.   This  is 
another  manifestation  of  the  negative  bias  in  the  prediction 
of  pressure  forces  in  the  trough  region  of  waves  which  is 
introduced  by  the  stretching  of  the  vertical  coordinate 
dependent  function.   Neither  the  linear  nor  the  nonlinear 
simulations  resulted  in  pressure  force  sequences  with  kurtosis 
measures  greater  than  0.50.   The  low  measures  of  both  the 
skewness  and  kurtosis  indicate  that,  statistically,  the  simu- 
lated pressure  force  sequences  are  almost  Gaussian  in  con- 
trast to  the  measured  values  which  generally  demonstrate  low 
skewness  measures  but  significantly  higher  than  Gaussian  kur- 
tosis measures. 
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Finally,  the  changes  in  the  standard  deviation  and  kur- 
tosis  resulting  from  the  transformations  from  the  linear  to 
nonlinear  realizations  were,  with  one  exception,  negative 
while  the  changes  in  the  skewness  were  consistently  positive. 


CHAPTER  4 

CONCLUSIONS  AND  RECOMMENDATIONS 


1.   Random  Sea  Simulations  Correct  to 
Second  Perturbation  Order 


The  nonlinear  boundary  value  problem  which  describes  the 
propagation  of  irregular  random  surface  gravity  waves  in  water 
of  finite  depth  has  been  solved  correct  to  second  perturbation 
order.   The  expansion  of  the  random  velocity  potential  in  a 
perturbation  power  series  was  assumed  to  be  valid  for  both 
continuous  and  discrete  spectra.   The  linear  first  order  solu- 
tion was  determined  from  a  well  posed  Sturm-Liouville  problem 
in  the  vertical  coordinate  which  resulted  from  the  separation 
of  variables  in  a  rectangular  cartesian  coordinate  system. 
The  homogeneous  Neumann-type  bottom  boundary  condition  (BBC) 
and  homogeneous  mixed  boundary  condition  representing  the  com- 
bined free  surface  boundary  condition  (CFSBC)  resulted  in  a 
deterministic  eigenvalue  problem  from  the  linear  spectral 
operator  which  is  unique  to  within  an  arbitrary  multiplicative 
constant. 

The  second  order  solution  was  found  to  be  a  function  of 
the  linear  solution  through  the  combined  free  surface  boundary 
condition.   The  separation  constant  for  the  second  order  solu- 
tion was  found  to  be  the  sums  and  differences  of  the  linear 
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eigenvalue    spectrum.      The    spectral    operator   for  the    second 
order  solution  was    computed   and  was    found  to  be  nonsingular. 
The    trivariance    function  was    shown   to  be   the   best  measure    of 
the   second   order   correction    to   the    linear  solution.      The 
coarse   measure    of   the    trivariance    function,    the    skewness,    was 
used  to   determine    the   measure   of  the    second  order   corrections 
as   well    as    the   non-Gaussian  nature    of   the    linear   simulation. 
Recognizing   that   the    second  statistical   moment   is   not    closed 
to   the   second  perturbation   order,    the    second   order   spectrum 
was    computed   and   found   to  be    the    convolution   of   first   order 
Gaussian    spectral   densities.      This    important   result    for   a 
linear  Gaussian   first    order   spectrum  has   been    given  by   Tick 
[123]    and   Longuet-Higgins    [79]. 

Random  realizations    of   four  measured  hurricane  wave 
records   were    simulated   correct    to   second  order.      Each    of   the 
four  measured   records   was    simulated  using   the    smoothed  mea- 
sured  spectrum  and  the    Bretschneider   spectrum  with    an  equiva- 
lent  variance    and  best    least-squares    fit   to   the   peak    fre- 
quency.     Three    linear   and   three   nonlinear   random  realizations 
for  each  measured  hurricane    record  were    simulated   from  both 
the    smoothed  measured  spectrum  and  the    Bretschneider   spectrum. 
The    root-mean-square   error  between   the   measured   and  Bret- 
schneider  spectra  were    lowest   for  three    relatively  broad 
spectra   and  was    the   highest   for   the    one    relatively  narrow 
spectrum. 

Each   of  the    linear   realizations   was      found   to  be   slightly 
non-Gaussian.      Each   of   the   second  perturbation   order 
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realizations  made  a  positive  contribution  to  the  linear  skew- 
ness  measures.   The  cumulative  probability  distributions  from 
the  simulations  demonstrated  good  agreement  with  both  the 
measured  distributions  and  a  Gaussian  distribution  with  the 
same  mean  and  variance  within  two  standard  deviations.   Out- 
side of  two  standard  deviations,  the  agreement  varied  consid- 
erably.  The  largest  value  of  the  nonlinear  realizations  were 
approximately  0.5  to  2.0  standard  deviations  greater  than  the 
corresponding  largest  values  from  the  linear  simulations. 
Distributions  of  realizations  synthesized  from  the  smoothed 
measured  spectra  demonstrated  the  same  variability  outside 
two  standard  deviations  as  distributions  of  realizations  syn- 
thesized from  the  Bretschneider  spectra.   This  indicates  the 
importance  of  the  random  phase  angles  in  the  distributions  of 
random  surface  gravity  waves  and  the  importance  of  outlier 
spectral  components  in  nonlinear  interactions. 

In  general,  the  distributions  of  nonlinear  realizations 
from  the  Bretschneider  spectrum  were  as  good  as  those  from  the 
smoothed  measured  spectrum  and  indicates  that  it  is  valid  for 
design  use  in  simulating  nonlinear  hurricane  waves  correct  to 
second  order. 


2 .   Wave  Induced  Pressure  Forces  on  a 
Vertical  Piling  Computed  by  Digital  Linear  Filter 
with  "Stretched"  Vertical  Coordinate 


Three  different  types  of  random  realizations  were  used 
as  inputs  to  the  digital  linear  filter  technique  to  predict 
the  kinematics  required  to  estimate  pressure  forces  using  the 
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Morison  equation.      The    three    types   of   realizations   were: 
(1)    the   measured  sea   surface    realization;    (2)    the    linear   and 
nonlinear   realizations    simulated   from  the   smoothed  measured 
spectra;    and    (3)    the    linear   and  nonlinear  realizations    simu- 
lated  from  the    Bretschneider   spectra  having   an   equivalent 
variance    and  best    least-squares    fit   to   the    smoothed  measured 
spectra.      None   of  the   values    of  the    root-mean-square   pressure 
from  the    three   predicted  pressure    force    realizations   were 
equal    to   or  exceed   the   measured  pressure    force   values.      The 
predicted  pressure    force    sequences    computed   from  the   measured 
sea   surface    realizations   using   the   modified  Dean   and  Aagaard 
force    coefficients    always   yielded   the   highest   predicted   root- 
mean-square   pressures    and  demonstrated   a   spectral    shape   which 
approximated  very   closely   that    of   the   measured  pressure    force 
spectra.      This    indicates    that   the    digital    linear   filter   tech- 
nique   is    appropriate    for   analysis    and  design    for  waves   with 
amplitudes    that    are    at    least   one -half   the   breaking   amplitude 
provided   the    appropriate    random  realization   is    filtered. 

The    realizations    simulated   from  the    Bretschneider   spec- 
trum using   only   one    seed  number   for   the    random  number   genera- 
tor  consistently  yielded  root-mean-square  pressure    forces 
which  were   slightly  higher  than   those   realized   from  the 
smoothed  measured  spectra.      The    skewness    of  the   pressure 
force    sequence    computed   from  the    linear   realizations    from 
both   the    smoothed  measured   and  Bretschneider   spectra  were 
consistently  negative.      This   negative   bias    is    the    result   of 
the   vertical    coordinate    stretching   function.      The    skewness 
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measure  of  the  pressure  force  sequence  computed  from  the  non- 
linear realizations  from  both  spectra  were  consistently 
slightly  positive.   This  indicates  that  the  positive  skewness 
of  nonlinear  sea  surface  realizations  tends  to  compensate  for 
the  negative  skewness  introduced  by  the  vertical  coordinate 
stretching  function.   Pressure  force  simulations  using  the 
digital  linear  filter  in  conjunction  with  the  vertical  coordi- 
nate stretching  function  should  be  obtained  from  nonlinear 
realizations  with  a  positive  skewness  measure  in  order  to 
ensure  that  the  pressure  force  sequences  are  not  negatively 
biased. 

Three  of  the  measured  pressure  force  sequences  demon- 
strated slightly  positive  skewness  measures.   The  one  sequence 
with  a  slightly  negative  pressure  force  skewness  measure  was 
obtained  from  a  record  with  a  slightly  negatively  skewed  sea 
surface  realization.   Although  the  skewness  measure  should 
theoretically  be  positive  for  slightly  non-Gaussian  processes, 
the  sign  of  the  skewness  of  the  measured  pressure  force 
sequence  corresponded  to  the  sign  of  skewness  of  the  measured 
sea  surface  realization. 

The  kurtosis  of  the  measured  pressure  force  sequences 
from  the  relatively  narrow  band  Record  No.  06886/2  was  found 
to  be  platykurtic  (i.e.,  K  <  3) ,  while  the  three  measured 
pressure  force  time  sequences  were  found  to  be  leptokurtic 
(i.e.,  K  >  3).   The  measured  pressure  force  sequences  were 
found  to  be  nearly  Gaussian  between  two  standard  deviations 
as  were  the  sea  surface  simulations.   The  simulated  pressure 
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forces  synthesized  from  the  two  types  of  spectra  were  found  to 
agree  with  the  measured  pressure  force  realizations  between 
two  standard  deviations.   For  values  outside  two  standard 
deviations,  the  agreement  varied  considerably. 

3.   Recommendations  for  Additional  Applications 

Chu  and  Mei  [35,36]  have  shown  that  the  application  of 
the  method  of  multiple  scales  (cf.  Neyfeh  [95])  to  slowly 
varying  Stokes  waves  introduces  a  correction  to  the  dispersion 
equation  (and,  therefore,  to  the  eigenvalues)  at  second  order 
in  the  perturbation  expansions.   The  extension  of  the  non- 
linearization  procedure  for  simulating  nonlinear  random  waves 
by  the  application  of  the  method  of  multiple  scales  to  a 
spectrum  of  random  waves  would  offer  the  advantage  of  intro- 
ducing a  nonlinear  correction  to  the  second  order  interaction 
matrix  through  the  second  order  corrections  to  the  eigen- 
values.  The  extension  would  be  more  complicated  since  the 
nonlinear  boundary  value  problem  would  have  to  be  solved  cor- 
rect to  third  order  (cf.  Chu  and  Mei  [35,36]).   In  addition, 
the  spectral  amplitudes  and  phases  would  become  slowly  varying 
functions  of  time  and,  consequently,  nonstationary.   Nonsta- 
tionary  spectral  analyses  involve  more  complex  evolutionary 
spectral  techniques  (cf.  Brown  [31]). 

Hasselmann  [60]  has  developed  a  general  solution  for  the 
nonlinear  interaction  of  shear  currents  and  surface  gravity 
waves.   This  general  solution  assumes  that  the  fluid  is 
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inviscid  and  that  the  velocity  field  may  be  decomposed  into 
rotational  and  irrotational  components.   The  vertical  vor- 
ticity  distribution  is  expanded  in  a  complete  orthonormal  set 
over  the  finite  depth,  h,  given,  e.g.,  by  the  following  series 

^r  ="Vf  ^  A^-f(t)-sin{^  (h+z)}exp  i{^-x}  (3.1) 

n 

The  modifications  of  the  Bernoulli  equation  due  to thepresence 
of  the  shear  current  are  given  by  Hasselmann  in  Appendix  C 
to  [60].   In  addition,  the  Eulerian  equations  for  nonlinear 
gravity  waves -shear  current  perturbations  are  reduced  to  the 
general  Hamiltonian  equations  by  Hasselmann  [60]  and  some  of 
the  difficulties  encountered  in  applications  are  discussed. 
The  nonlinearization  procedure  developed  in  this  study  could 
be  extended  to  include  shear  currents  by  this  method  but  the 
algebraic  details  may  prove  to  be  quite  complex. 

The  stretching  function  utilized  to  evaluate  the  hori- 
zontal kinematic  fields  introduced  a  negative  bias  compared 
to  values  computed  from  the  unmodified  linear  theory  velocity 
potential.   Alternate  stretching  functions  should  be  investi- 
gated and  compared  with  the  kinematic  field  values  predicted 
for  vertical  elevations  above  the  mean  sea  level  computed  by 
a  Taylor  series  expansion  similar  to  the  method  described  by 
Reid  [106]. 

The  "peak  enhancement  factor"  (cf.  Hasselmann  et  al. 
[62])  was  not  employed  for  the  random  simulations  synthesized 
from  the  Bretschneider  spectrum.   The  application  of  this 
factor  would  be  of  value  in  realizing  larger  nonlinear  self- 
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interaction   contributions    from  the   enhanced  spectral    ampli- 
tudes  near   the   peak   frequency   of  the    spectrum   compared   to   the 
nonlinear  self-interaction   contributions    from  a   Bretschneider 
spectrum  having  equal   variance.      The    application   of   the   peak 
enhancement    factor   results    in   transforming   the    Bretschneider 
spectrum   from  a   two-parameter   spectrum  to   a   five-parameter 
spectrum.      The    correlation   of  these    five   parameters   with 
parameters   measurable    from  measured  spectra  has   not  been   com- 
pleted   (cf.    Hasselmann   et^  al.     [62]).      The   transformation   to   a 
five -parameter   spectral    representation   still   results    in   a 
smoothly   continuous    spectrum   from  which   outlier   spectral   esti- 
mates which    are    capable    of  large    self-interaction   contribu- 
tions  have  been   effaced. 

Finally,    Green's    functions   have   been   extremely   useful    in 
determining  pressure    forces    on    fixed   and   floating   obstacles 
for   the    linear   deterministic  boundary  value   problem    (cf. 
Hildebrand    [65]    and  Wehausen    [127]).      Free    space    Green's    func- 
tions  have    also  been   useful    in   evaluating  wave   propagation 
problems    for   random  boundary  value   problems    (cf.    Boyce    [26] 
and  Frisch    [51]).      The    random  linear  velocity  potential   may 
be   expressed   as    a   spatial    integral    over   the   boundaries   with 
the   Green's    function   as    the    resolvent   kernel    of  the    integral. 
If  this    linear   integral    solution   could  be   successfully   em- 
ployed  in   the    inhomogeneous    forcing   terms   of  the   second  order 
combined   free   surface   boundary   condition,    the    result  would  be 
an  extension   of  the    Green's    function   to   second  order   and  would 
be   extremely   useful    in   radiation,    scattering   and   diffraction 
problems. 


APPENDIX  A 
THE    FOUR- PRODUCT   MOMENT    FOR  A   GAUSSIAN  VARIATE 

The    factoring   of   the    fourth   order   statistical   moment   of 
a  jointly  normal   Gaussian  variate  was   used  by   Tick    [123]    and 
Longuet-Higgins    [79]    to   obtain   irreducible   products   of  co- 
variances.      Parzen    [98;    p.    93]    and  Middleton    [88;    p.    343] 
give    different  methods    for   obtaining   this    important    and  useful 
result.      The    derivation   given  below   is    included   for   conven- 
ience. 

The    four-product   moment    of  the    four   random  variables 
Z- ,    Z^,    Z_,    and   1.    is    to  be   evaluated   from  the    following 
expression: 


'{l^Z^l^l^}     =    E{l^l^}'E{l^l^}+E{Z^l^'\-E{l2l^)  + 


*e{Z^1^}-e[1^1^} 


(A.l) 


Consider   the    random  variable    Z    as    a   four-element    column 
vector   given  by 

'Z. 

J. 

Z  =   .2  (A.2) 

with  vector  matrix  transpose  given  by 


[Zi,Z2,Z3,z,; 


(A.  3) 
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and  with  zero  mean 
e{Z}    =   0 


(A.4) 


The  characteristic  function,  <i>(u)  ,    is  the  Fourier-Stieltjes 
transform  of  the  distribution  function  of  the  random  variable, 
i.e.  , 

00 

(j)„(u)  =  /  exp  i{u  z}dP(z) 

^        -co 

<t)2(u)  =  E{exp  i{u'^Z}} 

1   T 

(t)^  (u)  =  expf-j  u  -P^'u} 


where  Py  is  the  covariance  matrix 


P^  =  e{ZZ'} 


e{Z^Z^}  £{1^1^}  EiZ^Z^}  ^{Z^Z^} 

e{Z^Z2}  £{^2^2^  ^^^3^2^  e{Z4Z2} 

EiZ^Z^}  EiZ^Z^}  e{Z^Z^}  EiZ^Z^} 

^{Z^Z^}  e{Z2Z4}  eCZ^Z^}  ^{Z^Z^} 


(A.  5) 


(A.  6) 


Upon  carrying  out  the  matrix  multiplication,  the  following 
result  is  obtained: 

4   4 


)(u,  ,u-,u_,u.)  =  exp{-y  E    E   u.u.  b{Z.Z.} 
z.-|.-|i3     ij 

17       7       7  J 

^1»^2»^3'^4 


(A.  7) 


The  four-product  moment  may  be  computed  by  differenti- 
ating the  characteristic  function  with  respect  to  its  four 
arguments  and  evaluating  the  result  with  the  argument  set 
equal  to  zero,  i.e.,  ' 


Eil^l^l^l^-^    -    s      5^^u,9u 


(t)(0, 0,0,0) 


1      2""3""4        1^,1^,1^,1^ 
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(A.  8) 


Note    that,    in   general 


34.(u^=0) 


=    -    E      u.    £{Z    Z.} 


(A. 9a) 


so    that 


3(f)(0,U2,U2,u^) 


E      Uj^-f{Z^-Z2-(})(0,U2,U3,U4) 

7        7        7        7 


3    (f)(0,0,U3,U4) 

^l'^2'^3'^4 
3u,3u2 


(A. 9b) 


4  4 

=    [-E{Z^Z2}+    Z      u^E{Z2ZjJ-    E      Uj^   e{1^\}\ 


.(J)(0,0,U3,U4)  (A. 9c) 

^l'^2'^3'^4 


3-^({)(0,0,0,U4) 

^l'^2>^3>^4 
3u, 3U28U2 


u^£{Z^Z2}-£{Z3Z^}+U4e{Z2Z3}£{Z^Z4} 
_+U4e{Z2Z4}-£{Z^Z3}-u^b{Z2Z4}-e{Z^Z4}-£{Z3Z4} 


3  ^(|)  (0,0,  0,0) 

^l'^2'^3'^4 
3u,3u23u23u^ 


.(^(0,0,0,U4)  (A.9d) 

Zi,Z2,Z3.Z4 


=   e{Z^Z2}-e{Z3Z4}+e{Z2Z3}-e{Z^Z4}+ 


+e{Z2Z4}-£{Z^Z3}  (A.9e) 


which  is  the  desired  formula.  The  vanishing  o£  the  three- 
product  moment  for  a  Gaussian  variate  may  be  obtained  from 
Eq.     (A.9d)    by   setting   u^eO. 


APPENDIX  B 
ALTERNATE  QUADRATIC  FILTER 

The  realizations  correct  to  second  perturbation  order 
shown  in  the  figures  in  Chapter  3  were  synthesized  from  a 
convolution  and  linear  siommation  in  the  frequency  domain  in 
order  to  exploit  the  spectral  representation  of  differential 
operators  of  the  perturbed  set  of  linear  problems  by  means  of 
the  fast  Fourier  transform. 

The  general  applications  of  nonlinear  filtering  methods 
in  the  time  domain  for  frequency  modulated  problems  have  been 
given  by  Wiener  [132].   A  discussion  relating  to  the  general 
quadratic  filtering  problem  was  given  by  Tick  [124].   The 
application  to  the  simulation  of  nonlinear  random  surface 
gravity  waves  by  a  quadratic  filter  which  is  derived  in  this 
appendix  was  suggested  by  Borgman  [25]. 

The  time  series  realization  for  the  sea  surface  may  be 
expressed  by  the  following: 

00  CO 

n(t)   =  hQ  +  /   h^(T)-?(t-T)dT  +  //  h2(T^,T2)- 

-^Ti(t-T^)  •^n(t-T2)dT^dT2  (B.l) 

where  ^(t-x)  is  a  Gaussian  white  noise  process  such  that 
C(t)  ^   N(0,1)  (B.2) 
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and  ^ri  is  the  first  order  realization  given  by  the  first 
integral.  The  value  of  the  zeroeth  order  kernel,  hr,,  is 
determined  by  requiring  that  the  process  have  zero  mean,  i.e., 

OO  00* 

E{n(t)}     =     h^     +    /     i:i(T)-E{cCt-T)}dT     +    //     h2(T^,X2)' 

•E{^n(t-T-|^)  •-,^n(t-T2)}dTj^dT2 

=   0  (B.3) 

which   reduces   to 

OO  00 

hQ    =    -   /    dx,    /   h^Cx,  ,T2) •Y(T,-T2)dT2  (B.4) 

-co  -OO 

1^1^ 

by   virtue   of  Eq.    (B.2) . 

The    first   order  kernel    is    the    impulse    response    function 

of  the    deterministic   spectral    representation   given  by 

M 
h^Cx)    =    27T      Z        C   5(x-mx    )  (B.5) 

m=-M 

where    the   Fourier   coefficients    are    determined   from 

Cj^  =   ^   /     /SCajAa    -exp    i{moT^}da  (B.6) 

"^     0       nn 

and  a      is    the    cut-off   frequency. 

The    second   order  kernel    requires    a   double   Dirac   delta 

comb    and   is    given  by 

^      M  M 

h^ir^,^^)    =    4Tr        ^   ^        C        6(xi-nx^).6(x2-mx„)  (B.7) 

n=-M  m=-M        ' 

where  the  Fourier  coefficients  are  determined  by 
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^n,in  "   ~T    -^     "^    Hia^,a^)-exY>   i{na^T^+ma2TQ}da^da2  (B.8) 


Substituting   the    appropriate    impulse    response    functions,    the 

random  sea   realization   correct   to   second  perturbation   order 

becomes : 

M  °° 

n(t)    =    2tt      E        C^-  S    ?(t-T) -6  [x-mx    ]dT+ 
m=-M  -0°  ° 

^      M  M  oo 

n=-M  m=-M        '      -°° 
.6[T2-mT^]    dT^dT2  (B.9) 

The  convolution  given  by  Eq.  (B.9)  may  now  be  computed  in 
the  time  domain  from  the  appropriate  known  Fourier  series 
expansions  of  the  energy  density  spectrum  and  the  nonlinear 
interaction  kernel  derived  in  Section  4,  Chapter  2.   The 
linear  first  order  realization  given  by  the  first  integral  in 
Eq.  (B,9)  must  be  evaluated  first  and  then  the  convolution  of 
this  realization  with  the  Fourier  coefficients  representing 
nonlinear  interaction  matrix  is  computed  and  the  two  time 
series  added  in  accordance  with  Eq.  (B.9). 


APPENDIX  C 

EFFECT  OF  THE  HORIZONTAL  SPATIAL  SEPARATION 
BETWEEN  WAVE  STAFF  AND  INSTRUMENTED  PILING 
ON  THE  EVALUATION  OF  PRESSURE  FORCE  COEFFICIENTS 


The  effect  of  the  horizontal  spatial  separation  between 
a  wave  staff  and  an  instrumented  platform  piling  on  the  eval- 
uation of  pressure  force  coefficients  by  a  regression  of  the 
measured  and  predicted  values  can  be  most  significant. 
Figure  C.l  indicates  the  geometry  for  the  wave  staff  and 
instrumented  piling  for  the  Wave  Project  II  data  acquisition 
system.   The  coordinate  axes  are  aligned  with  the  platform. 
The  measured  force  may  be  expressed  by  the  Morison  equation 
as 

F^  =  Kp-u|ul.Kj-u^  (C.l) 

where  K^  is  a  constant  representing  the  drag  coefficient 
sea  water  density,  and  piling  geometry  and  K^  is  an  inertia 
constant  representing  the  inertia  coefficient  and  sea  water 
density,  and  piling  geometry.   The  predicted  pressure  force 
may  also  be  given  by  the  Morison  equation  as  follows: 

F^  =  K-u'lu'hK^ul  CC.2) 

Assuming  that  kinematics  from  linear  wave  theory  are  approxi- 
mately correct,  the  horizontal  velocity  and  acceleration  for 
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Direction  to 
Wove  Staff 


Wove     Staff  •(xy) 


Direction     of 
Wove     Propagation 


///////////////////////////X     1^ 


Instrumented 
Piling 


Waves 


Figure  C.l. 


Orientation  o£  Wave  Staff  and  Instrumented 
Piling  for  Wave  Project  II. 
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the  measured  force  are 

u  =  /A^  •  cos{at}  (C.3a) 

u^  =  -y/A^    •  sin{at}  •               (C.3b) 

and  for  the  predicted  force  are 

u'  =  /A^  •  cos{at-(k^-x^+ky.y^)}  (C.4a) 

u|  =  -/a^  •  sin{at-(k^-x^+ky-y^)}  (C.4b) 

where  '^^,1^    are  the  coordinates  of  the  wave  staff,  k  ,k   are 
the  wave  numbers  referenced  to  the  platform  coordinates,  and 
/A^,/a^  are  functions  of  the  vertical  coordinate  along  the 
pile  and  the  wave  amplitude. 

The  spatial  phase  shift  may  be  related  to  the  angular 
frequency  through  the  orthogonal  components  of  the  wave 
celerity  referenced  to  the  platform  coordinate  system  to 
yield 

k  -x   +  k  .y  =   o    '    {-^  ^   -^^    -   o    '    V  (C.5) 

X  s    y   s         X   ^y 

The  square  of  the  error  between  the  measured  and  pre- 
dicted forces  during  the  passage  of  a  wave  crest  may  be 
expressed  by 

9    t/4        7 
e^  =  /   {F„-F^}^  dt  (C.6) 

-T/4   "^  P 

The  presence  of  the  absolute  value  of  the  horizontal  velocity 
in  the  drag  force  term  requires  that  the  integral  be  decom- 
posed into  two  intervals  according  to  the  value  of  the  sign 
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of  the  trignometric  term  containing  the  phase  shifted  argu- 
ment.  The  cosine  becomes  negative  when 


t-f  <  -T/4    for  t  >  -T/4 
and  becomes  positive  when 
t-f  >  -t/4    for  t  <  t/4 


(C.7) 


CC.8) 


Therefore,  the  integral  in  Eq.  (C.6)  may  be  decomposed  into 

the  following  two  integrals: 

-T/4+t ' 
e^  =   /    {-Kt.-A,  cos'^(at)-Ky/A^    sin(at)- 
.T/4      ^   ^  ^   ^ 

-[-K^A^  cos^{a(t-t')}-KjA^  sin{a  (t-t ' ) }  ]  }^  dt+ 

+  /     ^K^-A^  cos^(at)-K./A^  sinCat)- 
-T/4+t'    ^   ^  ^      ^ 

-[^A^    cos^{a(t-t')}-K'/A^  sin{a(t-t)}]}^  dt        (C.9) 


Equation  (C.9)  is  to  be  minimized  with  respect  to  the  unknown 

force  constants,  K^  and  KJ.   Differentiating  Eq.  (C.9)  first 

with  respect  to  the  drag  coefficient,  K'   and  then  with 

respect  to  the  inertia  coefficient,  K'   and  equating  each 

result  to  zero  yields  the  following  two  equations: 

-T/4+t'  .  . 

/    [-Kp-A^-cos^(at)-Kj-/A^-sin(at)+K^-A^-cos^{a(t-t')}+ 

+KJ.-/A^-sin{a(t-t')}]-cos^{a(t-f)}dt+ 


t/4 
/ 
■T/4+t' 


Kp-A^-cos  (at)-Kj-/A^-sin(at) 
_-K^-A^-cos^{a(t-t')}+KJ-/A^-sin{a(t-t')} 


•cos^{0(t-f)}dt  =  0 


(C.lOa) 
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■T/4  +  t 

/ 
•T/4 


Kp-A^-cos  (at)-Kj.-/;i^-sin(at) 
.+K^-A^-cos^{a(t-t')}+KJ-/A^-sin{a(t-t')}_ 


sin{a(t-t')}dt 


T/4 
+  / 
-T/4+t' 


Kjj-A^-cos  (at)-Kj-/A^-sinCcrt) 
.-K^-A^-cos^{aCt-t')}  +  KJ-/;i^-sin{a(t-t' )} 


sin{a(t-t')}dt  =  0 


(C.lOb) 


Define 

the  following  ratios: 

Pi 

S-^z 

4 

(C.lla) 


(C.llb) 


(C.llc) 


After  very  tedious  algebra,  the  following  dimensionless  ratio 
between  the  computed  drag  coefficient,  C^,  and  the  actual 
drag  coefficient  may  be  determined: 


'D  _  4    .77 

C^  ~  377  ^2 


I  sin{2af}-at'-[^  -  J]cos  (2at ' ) }  + 


l^jl    32 
^  TP^-  f^sin(af) 


(C.12) 


and  the  dimensionless  ratio  between  the  computed  inertia 
coefficient,  CI,  and  the  actual  inertia  coefficient  is 

C 
I  i"I 


|p  I 
^  =  cos(at')-  -r-p^-p  •  ^  {sin(at')+2sin(2at')} 


(C.13) 


For  very  small  separation  distances  compared  to  wave  period, 
the  following  asymptotic  values  are  easily  found: 
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From  Fig.    C.l,    the    following   identity   is    obtained: 


at'    =  k^-x^   .  ky-y^ 

=  k'cos   3'r'cos   a  +  k'sin   S-r-sin  a 

=  kr   cos(6-a)  (C.15) 

Figures  C.2  and  C.3  give  values  of  dimensionless  drag  and 
inertia  coefficient  ratios  as  a  function  of  at'  in  terms  of 
the  following  dimensionless  distance  for  various  ratios  of 
drag  and  inertia  pressure  force: 

^  =  2_u  .^.cos(6-a)  (C.16) 

where  L  is  the  wave  length  determined  from  linear  theory. 
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Figure  C.2.   Effect  of  Dimensionless  Wave  Staff- Instrumented 
Piling  Separation  Distance  on  Drag  Coefficients 
Determined  from  Morison  Equation  with  Linear 
Theory  Kinematics  (L  =  linear  theory  wavelength) 


1S5 


6 
5 

^ '           '         y^' 

■    1           II          1          1 

4 

_                                ////'~~\  \\\ 

r'  =  27rr    cos  (^-a) 

3 

_                              /////            \\v^ 

- 

.2 
1     2 

1      ' 
o 

IV//^^^^^ 

i^LA 

1 " 

Iv^T^^— ^^^^^^ 

.5    -1 
1.3 

^^ 

\j 

1-4 
5 

-  \-~^ 

%f  -  ■" : 

-5 
-6 

1           1           1          1          1 

1        1      1       1       1   . 

•3.14  -2.61  -2.04  -1.57  -1.04  -0.52 


0.52      1.04  1.57   2.04    2.61     3.14 


Dimensionless     Separation      Distance 


(Radians) 


Figure  C.3.   Effect  of  Dimensionless  Wave  Staff- Instrumented 
Piling  Separation  Distance  on  Inertia  Coeffi- 
cients Determined  from  Morison  Equation  with 
Linear  Theory  Kinematics  (L  =  linear  theory 
wavelength) . 
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